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PREFACE 


/'The  main  objective  of  this  study  is  to  develop  a  numerical  model  for 
routing  water  and  sediment  on  small  agricultural  catchments.  Part  2  of 
this  report  presents  a  detailed  description  of  the  models — The  model/ is 
developed  on  a  general  basis  so  that  it  may  be  applied  to  any  agricultural 
catchment,  and  it  can  be  used  to  simulate  the  effect  of  different  land  uses 
on  the  water  and  sediment  yields  from  the  modeled  catchment.  ^Paxt 

C  h  f  *  < 

presents  results  from  validations  of  the  model  using  several  data  sets 
including  data  from  natural  catchments.  Input  data  and  computer  coding 
details  are  given  in  the  Addendum  to  this  report. 
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INTRODUCTION 


Alluvial  streams  are  dynamic  systems  that  continuously  change  their 
configuration  and  state  in  response  to  either  changes  in  the  natural 
environment,  or  perturbations  introduced  by  man's  activity.  Frequently, 
these  changes  conduce  to  alteration  of  the  stream-channel  stability,  which 
often  results  in  channel  migration  and  shoaling. 

Among  the  leading  causes  of  channel  instability  are  several  that  are 
intimately  associated  with  land -management  and  conservation  practices 
carried  out  on  the  upland  areas.  They  are  (a)  clearing  of  land  that 
removes  the  soil-protective  and  flow- retardant  ground  cover,  which  in  turn 
leads  to  increased  erosion  and  flood  peaks;  (b)  installation  of  reservoirs 
for  flood  protection  and  irrigation  control,  which  upset  the  water-sediment 
equilibrium  downstream  of  those  structures;  and  (c)  excessive  soil  erosion 
resulting  from  uncontrolled  sources.  The  combined  effect  is  an  aggregate 
flow  of  water  and  sediment  coming  from  a  variety  of  point  and  non-point 
sources  within  the  upstream  catchments.  This  aggregate  yield  acts  as  a 
time  and  space  dependent  loading  on  the  streams  draining  the  catchments. 
If  this  loading  becomes  quite  different  from  that  which  the  streams  have 
adjusted  to,  the  result  is  a  breakdown  in  the  stability  of  the  channel 
system. 

The  catchments  contributing  to  the  loading  of  any  given  channel  system 
exhibit  in  general  a  great  variety  of  soils,  vegetation,  and  land  uses.  In 
order  to  effectively  assess  the  impact  of  these  catchments  on  the  loading 
of  the  channel  system,  it  is  necessary  to  develop  improved  methods  for 
predicting  the  effects  of  alternative  land  managements  of  those  catchments. 
There  is,  therefore,  a  need  for  the  development  of  mathematical  models  so 
that  the  hydrology  of  a  catchment  can  be  simulated  and  the  effects  of 
various  management  practices  understood  and  predicted.  In  response  to  this 
need,  the  goal  of  the  present  study  is  the  development  of  a  prediction 
model  for  estimating  sediment  yield  from  agricultural  catchments.  In 
developing  the  model,  the  following  specific  objectives  were  considered: 
(i)  estimate  the  amount  of  soil  loss  from  specified  soil-source  units  with 
homogeneous  characteristics;  (ii)  estimate  the  amount  of  water  and  sediment 
transported  out  of  the  catchments  through  the  principal  drainage  networks; 
and  (iii)  estimate  the  rate  of  channel  aggradation  and  degradation  along 
the  flow  system.  The  model  is  oriented  towards  the  needs  of  the  Corps  of 
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Engineers  for  better  means  of  assessing  the  impact  of  land-management 
practices  on  stream  channel  behavior. 

Many  hundreds  of  papers  have  been  written  concerning  studies  on 
various  aspects  of  hydrology.  For  this  reason  it  is  quite  impossible  to 
summarize  the  previous  work  that  has  led  to  the  current  understanding  of 
the  hydrologic  cycle.  Reference  can  be  made  to  some  of  the  existing 
comprehensive  hydrology  books  (i.e.,  Chow,  1964),  and  to  a  number  of 
American  Society  of  Agronomy  Monographs  (Luthin,  1957;  Van  Shilfgaarde, 
1974;  Hagan,  Haise  and  Edminster,  1967)  and  reports  (Pierre  et  al.,  1966; 
Neilsen,  Jackson,  Cary,  and  Evans,  1972)  that  discuss  current  knowledge  of 
soil-water-plant  system.  Excellent  reviews  of  the  progress  made  during 
recent  years  in  several  hydrologic  subjects  have  been  presented  by  Schaake 
(1975),  Amerman  et  al.  (1975),  Johnson  and  Meyer  (1975),  and  Nordin  (1975). 
Only  because  of  this  accumulated  knowledge  is  the  proposed  project  even 
feasible.  As  a  better  understanding  of  the  hydrologic  cycle  and  the  basic 
physical  laws  governing  it  have  evolved,  they  have  been  synthesized  into 
more  rational  and  physically  based  models.  Finally,  the  accessibility  to 
high-speed  digital  computers  has  made  possible  the  development  of  detailed 
comprehensive  hydrologic  models. 

Because  the  physical  processes  governing  catchment  behavior  are  very 
complicated,  many  past  studies  have  utilized  regression  models.  However, 
it  is  difficult  to  predict  the  response  of  a  catchment  to  different  land- 
management  activities  using  regressional  methods,  because  these  methods  are 
based  on  the  assumption  of  time  and  space  invariability.  This  assumption 
almost  always  fails  to  be  valid  in  the  case  of  natural  catchments. 

A  second  type  of  models  includes  lumped  parametric  simulation  methods, 
such  as  the  TVA  Continuous  Daily-Streamf low  Model  (TVA,  1972).  These 
models  simulate  the  response  of  a  given  catchment  by  adjusting  a  number  of 
coefficients,  with  little  physical  significance,  using  data  collected  under 
certain  environmental  conditions.  The  impossibility  of  relating  those 
coefficients  to  a  different  set  of  environmental  conditions,  seriously 
restrict  the  use  of  these  models  for  predicting  the  response  of  ungaged 
catchments . 

A  different  class  of  models  embodies  the  distributed  process 
simulation  methods.  These  techniques  use  mathematical  descriptions  of  the 
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basic  hydrologic  processes  being  modeled,  and  their  interaction.  In 
addition,  this  approach  tends  to  minimize  the  number  of  adjustable 
parameters  and,  whenever  possible,  relates  them  to  physical  quantities  that 
can  be  readily  measured  in  the  field. 

The  Stanford  model  (Crawford  and  Linsley,  1966)  was  one  of  the  first 
general  models  developed  to  simulate  runoff  from  a  catchment.  It  is 
basically  a  lumped-parameter  model,  although  large,  heterogeneous 
catchments  can  be  subdivided  into  subcatchments  if  sufficient  data  are 
available  to  define  model  parameters.  The  model  has  gained  widespread  use 
and  as  a  result  has  undergone  numerous  modifications.  Holtan  and  Lopez 
(1970)  have  described  the  USDAHL-70  model  of  catchment  hydrology.  Although 
this  model  is  basically  lumped,  a  heterogeneous  catchment  can  be  broken 
down  into  smaller  homogeneous  areas.  An  attempt  is  made  to  incorporate 
spatial  variability  by  dividing  the  catchment  into  land  capabilities 
classes  that  correspond  to  uplands,  hillslopes,  and  bottom  lands.  Dawdy  et 
al.  (1972)  reported  on  a  lumped-system  model  similar  to  the  Stanford  model 
which  describes  surface  runoff  from  small  catchments.  TVA  (1972)  recently 
described  a  lumped  daily-streamf low  model  with  sixteen  parameters,  five  of 
which  require  optimization.  This  model  has  been  reasonably  successful  in 
predicting  daily  streamflows. 

A  continuous  distributed  model  is  not  yet  available.  However,  several 
single-event  distributed  models  that  include  part  of  the  hydrologic  cycle 
have  been  introduced  since  the  pioneering  works  of  Wooding  (1966)  and 
Woolhiser  and  Liggett  (1967).  Since  then,  a  cascade  of  various  sizes  and 
slopes  (Brakensiek,  1967;  Kibbler  and  Woolhiser,  1970)  or  converging 
inverted  cone-shaped  surfaces  (Woolhiser,  1969)  have  been  used  for 
geometric  representation  of  complex  topographies.  The  works  of  these  and 
other  investigators  have  let  to  the  acceptance  of  the  kinematic-wave 
approximation  as  an  adequate  model  of  shallow  overland  flow  and  flow  in 
channels.  The  reductionist  approach  to  watershed  simulation  was  introduced 
by  Huggins  and  Monke  (1970).  They  employed  a  square  grid  for  decomposing  a 
complex  catchment  into  elemental  surface  units.  Most  physically-based 
overland  flow  models  used  simplified  lumped-system  infiltration  models. 
Smith  and  Woolhiser  (1971)  were  the  first  to  introduce  a  distributed 
infiltration  model,  derived  from  soil  moisture  flow  theory,  to  calculate 
point  infiltration  rate,  and  therefore  rainfall  excess  rate.  The  foregoing 


concepts  have  been  incorporated,  in  one  way  or  another,  into  more  detailed 
models  recently  reported  by  Simons  et  al.  (1975)  and  Smith  (1976).  These 
models  use  flow  routing  techniques  based  on  the  kinematic-wave 
approximation  of  the  flow  governing  equations.  Since  its  formulation  by 
Lighthill  and  Whitham  (1955),  the  kinematic  wave  approximation  has  received 
extensive  application  to  catchment  runoff  modeling.  This  approximation  is 
restricted  by  the  assumption  that  the  friction  slope  equals  the  stream  bed 
slope,  but  it  has  been  found  to  be  applicable  in  many  stream  flows  and  in 
most  overland  flow  situations.  In  addition,  the  kinematic-wave  formulation 
admits  an  analytical  solution  by  the  method  of  characteristics  (Eagleson, 
1970;  Kibler  and  Woolhiser,  1970;  Li  et  al.,  1975b;  Singh,  1975).  This 
analytical  solution  has  two  main  advantages  over  other  numerical  solutions. 
It  eliminates  the  wave-celerity-damping  and  phase  lag  usually  induced  by 
numerical  schemes;  and,  in  addition,  results  in  faster  computational 
procedures.  In  spite  of  these  advantages,  applications  of  this  analytical 
solution  have  been  restricted  in  the  past  to  catchment  models  with  a  high 
degree  of  geometric  abstraction.  The  reason  is  the  formation  of  kinematic 
shock  waves  (Lighthill  and  Whitham,  1955;  Kibler  and  Woolhiser,  1970; 
Harley  et  al.,  1970;  Whitham,  1974).  Formally,  innumerable  shock  waves  can 
be  generated  during  the  routing  process,  as  a  result  of  the  time  and 
spatial  discretization  of  precipitation  and  the  physical  characteristics  of 
the  catchment.  In  the  past,  the  existence  of  these  shocks  has  frequently 
been  ignored  by  using  approximate  numerical  techniques.  This  practice, 
however,  may  not  be  considered  as  valid  particularly  when  the  foundation 
and  the  physical  relevance  of  the  kinematic  wave  approximation  is  under 
investigation.  It  is  well  known  that  shock  formation  is  intrinsic  to  the 
hyperbolic  equation  governing  kinematic  theory.  Further,  they  are 
considered  to  be  the  manifestations  of  higher  order  effects  such  as 
formation  of  monoclinal  flood  waves,  bores,  etc.  These  discontinuities 
play  important  roles  in  the  dynamics  of  hydraulic  systems  and  an  ad  hoc 
smoothing  by  numerical  means  does  not  necessarily  make  the  theory  look 
better.  The  model  described  in  the  present  report  introduces  a  new 
solution  to  the  kinematic  approximation,  which  retains  the  dynamic  effects 
of  the  shocks  by  routing  the  discontinuities  as  they  appear.  Certain 
simplifying  assumptions  are  made  which  permit  closed  form  solutions  and  an 
efficient  numerical  algoritlim,  based  on  the  method  of  characteristics.  The 
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resulting  procedure,  called  an  approximate  shock  fitting  scheme,  preserves 
the  effect  of  the  shocks  without  the  usual  computational  complications  and 
compares  favorably  with  existing  finite  difference  solutions  (Borah  et. 
al.,  1980). 

Different  types  of  sediment  production  models  have  appeared  widely 
dispersed  in  the  technical  literature.  Reference  can  be  made  to  a  recent 
review  presented  by  Heinemann  and  Piest  (1975)  and  to  a  publication  of  the 
Agricultural  Research  Service  (1975).  Several  regression  equations  for 
predicting  gross  soil  erosion  have  been  proposed.  The  most  commonly  used 
among  these  is  the  so-called  universal  soil  loss  equation  (USLE)  proposed 
by  Wischmeier  and  Smith  (1965).  Other  equations  of  similar  nature  have 
been  developed  by  Musgrave  (2947),  and  Gottschalk  and  Brune  (1950).  In 
these  equations,  the  soil  loss  rate  is  correlated  with  storm,  land,  and 
vegetation  characteristics.  Such  equations  are  applicable  on  seasonal 
basis  or  longer.  Also,  they  do  not  take  advantage  of  the  physical 
processes  occurring  within  the  catchments;  hence,  it  is  not  possible  to  use 
them  on  large,  complex  basins.  Williams  (1972)  modified  the  USLE  to  make 
it  applicable  for  predicting  storm  sediment  yields.  Onstad  and  Foster 
(1975)  combined  a  different  modification  of  the  USLE  with  the  USDAHL-70 
catchment  model  to  predict  sediment  yield  for  single  storms.  They  applied 
their  model  to  two  small  catchments  with  limited  success. 

The  first  physically-based  sediment  yield  model  was  reported  by  Negev 
(1967).  This  model  uses  the  Stanford  model  for  the  water  phase,  and  takes 
into  consideration  rainfall  soil  splash,  entrainment  by  overland  flow,  and 
rilling  and  gullying,  along  with  separate  channel  transport  of  fine  and 
coarse  sediment.  Sediment  production  is  evaluated  in  terms  of  power 
functions  of  water  discharge  containing  a  number  of  parameters  that  must  be 
calibrated.  A  modified  version  of  Negev's  model  has  been  incorporated  in 
the  Agricultural  Runoff  Management  (ARM)  model  recently  reported  by 
Donigian  and  Crawford  (1976).  The  aforementioned  models  of  Simons  et  al. 
(1975)  and  Smith  (1976)  also  incorporate  the  capability  of  describing 
sediment  movement  on  a  catchment  as  a  time  and  space  distributed  process. 
The  structures  of  these  two  models  are  similar;  however,  there  are 
differences  in  numerical  techniques  and  functional  relationships.  The 
sediment  movement  is  described  by  linking  the  excess-rainfall  flow 
equations  to  the  sediment  continuity  equation,  with  relations  describing 
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sediment  detachment  and  transport  capacity  at  any  point  on  the  surface  or 
in  a  channel.  A  similar  structure  has  been  incorporated  in  the  erosion  and 
sedimentation  component  of  the  present  model.  In  addition,  sediment  is 
routed  using  a  sediment  characteristic  scheme  that  takes  advantage  of  the 
efficient  analytical  solution  mentioned  above. 

Part  2  of  this  report  provides  a  detailed  description  of  the  model. 
Part  3  discusses  the  validation  of  the  model  on  sets  of  laboratory  and 
field  data.  Input  data  and  coding  details  are  given  in  Addendum  1.  This 
report  is  based  on  material  presented  in  an  earlier  study  by  Borah  (1979). 
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MODEL  FORMULATION 


The  model  basically  consists  of  two  intertwined  models:  one 
describing  the  hydrology  of  the  basin;  the  other  describing  the  associated 
erosion  and  sedimentation  processes.  It  simulates  the  movement  of  water 
and  sediment  as  a  time  and  space  distributed  process,  and  it  has  the 
ability  of  distinguishing  between  overland  and  channel  flows.  The 
catchment  is  regarded  as  consisting  of  a  mosaic  of  individual  subcatchments 
interconnected  by  channel  reaches.  The  model  can  thus  be  regarded  as  a 
cascading  process  in  which  the  output  of  one  or  more  subcatchments  becomes 
the  input  to  another  subcatchment  or  channel  reach.  In  mimicking  the 
overland  movement  of  water  and  sediment,  the  model  simulates  processes  of 
interception,  infiltration,  runoff,  detachment,  transport,  and  deposition 
of  sediment.  The  water  and  sediment  reaching  the  streams  are  then  routed 
through  the  channel  system,  and  the  rates  of  channel  aggradation  and 
degradation  are  computed.  The  basic  structure  of  the  model  is  graphically 
illustrated  in  Fig.  1.  The  details  of  the  model  components  are  given  in 
the  following  sections.  The  applicability  of  the  model  is  restricted  to 
catchments  where  the  streamflows  are  ephemeral,  the  subsurface  flow  and 
ground  water  movement  are  not  significant,  and  the  kinematic  wave 
approximation  for  flow  routing  is  valid.  The  computer  program  has  been 
written  assuming  a  uniform  distribution  of  rainfall.  However,  the  program 
can  be  easily  modified  to  accomodate  any  other  aerial  rainfall 
distribution. 

The  catchment  is  segmented  into  subcatchment  and  channel  reaches  to 
account  for  the  lack  of  uniformity  in  terrain,  soil,  and  land  use 
characteristics  in  most  natural  catchments.  Within  each  of  the  segments 
these  characteristics  are  treated  as  being  uniform.  The  subcatchments  are 
replaced  by  sloping  rectangular  areas  with  representative  length,  slope, 
width,  soil,  and  vegetative  characteristics.  The  channel  segments  are 
described  by  representative  cross-sectional  shape,  slope,  length,  and 
roughness.  Gravity  flow  logic  is  used  to  determine  the  computational 
sequence  as  explained  in  Appendix  I.  The  input  data  required  by  the  model 
includes  storm  characteristics,  geometry  data,  vegetative  cover  data,  soil 
data,  and  water  and  sediment  routing  data.  Details  on  input  data 
preparation  are  given  in  Addendum  1 . 
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Fig.  1.  Structure  of  the  model 
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2.1  INTERCEPTION.  NET  RAINFALL  RATE 

The  rainfall  excess  is  the  rainfall  contributing  to  the  water  flowing 
over  the  surface  of  an  overland  unit.  This  is  the  resultant  rainfall  after 
incurring  the  losses  due  to  interception,  evaporation,  transpiration  and 
infiltration.  At  the  beginning  of  a  rainfall  event,  some  rainfall  is  lost 
due  to  interception  at  the  canopy  and  ground  covers  and  thus  interception 
is  the  first  concern  in  computing  rainfall  excess.  The  interception 
component  adopted  in  this  model  is  based  upon  the  approach  proposed  by 
Simons,  et  al.  (1975). 

The  canopy  cover  and  the  ground  cover  are  the  two  major  features  which 
influence  the  motion  of  raindrops  before  reaching  the  ground  surface.  The 
canopy  cover  and  the  ground  cover  are  represented  by  the  canopy  cover 
density  and  the  ground  cover  density,  respectively.  The  canopy  cover 
density,  Dc,  is  defined  as  the  ratio  of  the  area  covered  by  trees  to  the 
total  area,  and  the  ground  cover  density,  D^,  is  the  ratio  of  the  ground 
area  covered  with  litter,  rock,  grass,  etc.,  to  the  total  area.  The  canopy 
cover  and  the  ground  cover  densities  are  assumed  uniform  on  a  flow  unit. 

Net  rainfall  is  the  quantity  of  rainfall  reaching  the  ground  after 
passing  through  those  covers.  Without  any  obstruction  like  canopy  cover  or 
ground  cover,  the  rain  reaches  the  ground  surface  at  the  same  intensity 
without  any  recognizable  loss.  In  the  presence  of  trees,  a  portion  of  the 
rainfall  is  stored  in  the  canopy  and  the  remainder  passes  through  the  trees 
as  throughfall  and  stemflow  (Zinke,  1965).  In  the  presence  of  ground 
cover,  a  portion  of  the  throughfall  is  intercepted  and  the  remainder 
finally  reaches  the  ground  as  net  rainfall. 

Let  I  be  the  rainfall  rate  (or  intensity)  at  time  t  at  a  level  above 
the  tree  canopy  as  shown  in  Fig.  2.  Let  I  be  the  rate  at  which  rain  is 
being  stored  in  the  canopy  at  time  t.  Then,  the  rainfall  rate  under  the 
tree  canopy  is  reduced  to  the  throughfall  rate  and  the  weighted  average 
throughfall  rate  is 

I  =  I  -  D  I  ,  (1) 
o  c  c 

where  stemflow  has  been  neglected.  Similarly,  let  I  be  the  rate  at  which 
rain  is  being  stored  in  the  ground  cover  at  time  t.  Then  the  weighted 
average  net  rainfall  rate  reaching  the  ground  is 

I  =  I  -  D  I  (2) 
n  o  g  g 
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2  Control  volumes  for  tree  canopy  and  ground  cover 


1 . 1-S 
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According  to  Horton  (1919),  the  total  interception  equals  leaf  storage 
capacity  plus  evaporation  loss  during  the  storm.  Zinke  (1965)  indicated 
that  .  .  usually  for  a  storm,  there  is  an  initial  period  during  which 
the  vegetation  cover  is  wetted  and  a  so-called  interception  storage 
capacity  is  satisfied.  This  is  followed  by  loss  from  this  storage,  and  the 
loss  is  dependent  upon  the  evaporation  opportunity  during  the  remainder  of 
the  storm." 

Let  be  the  interception  storage  capacity  of  a  tree  canopy  per  unit 
area  of  tree  canopy,  the  interception  storage  capacity  of  the  ground 
cover  per  unit  area  of  ground  cover,  E  the  mean  evaporation  rate  form  the 
interception  storages,  Sc  and  the  ratios  of  the  evaporating  surface  area 
to  the  horizontal  projected  area  for  a  tree  canopy  and  for  a  typical  ground 
cover,  respectively,  and  I  the  initial  interception  storage  content  which 
is  defined  as  the  ratio  of  the  initial  storage  to  the  interception  storage 
capacity.  Accordingly, 


Ic  =  I,  if  ,  /  K*)dt  S  (1  -  Is)  Vc  , 


I  =  ES  ,  if  ,  f  I(t)dt  >  (1  -  I  )  V 
c  c  J  v  s  c 

o 


Similarly, 


I  =  I  ,  if  ,  J  I  (T.)d  <  (1  -  I  )  V 
go’  to  s'  g 


I  =  E  S  ,  if  ,  J  I  (t)dt  >  (1  -  I  )  V 
g  g  o  °  S  8 


Mm*™?,..,  -"tfr  ' 


For  easy  handling  of  data,  a  parameter  rv  is  introduced,  defined  as 

the  ration  of  the  interception  storage  capacity  of  a  typical  canopy  cover 

to  that  of  the  ground  cover.  Therefore, 

V  =  r  V 
c  v  g 

and 

S  =  r  S 
c  v  g 

Eqns .  3  and  4  can  be  rewritten  in  the  discrete  form  as 

I  =  I  ,  if  I  IAt  <  (1-1  )r  V  ,  (7) 

c  S  V  g  ’  '  ' 

I  =  E  R  S  ,  if  I  IAt  >  (1-1  ) r  V  .  (8) 

c  v  g  s'  v  g  ' 

Similarly,  one  obtains 

I  =  I  ,  if  I  I  At  <  (1-1  )V  ,  (7’) 

go’  o  s'  g 

I  =  E  S  ,  if  I  I  At  >  (1-1  )V  .  (8’) 

8  g  o  s'  g 

Eqns.  1,  2,  7,  8,  V  and  8*  are  used  to  compute  the  discretized  net 

rainfall  rate. 

2.2  INFILTRATION.  RAINFALL  EXCESS  RATE 

Water  reaching  the  ground  surface  at  the  beginning  of  a  storm  passes 
through  into  the  soil  because  in  nearly  all  cases  the  initial  infiltration 
capacity  rate  is  greater  than  the  initial  net  rainfall  rate.  In  that  case, 
the  infiltration  rate  is  equal  to  the  net  rainfall  rate  and  there  is  no 
runoff.  Once  the  net  rainfall  rate  exceeds  the  infiltration  rate,  the 
excess  rainfall  either  runs  off  or  is  stored  on  the  ground  surface  as 
depression  storage  and  thus  ponding  takes  place.  The  time  elapsed  until 
the  beginning  of  excess  rainfall  is  known  as  the  ponding  time. 
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It  is  very  difficult  to  describe  mathematically  the  process  of  runoff 
originating  from  rainfall  excess,  because,  very  little  is  known  concerning 
the  magnitude  of  depression  storage.  Meaningful  observation  of  depression 
storages  are  not  easily  obtained.  Thus,  depression  storage  is  usually 
combined  with  interception  and  treated  as  an  initial  loss  with  respect  to 
storm  runoff  (Linsley  et  al.  1958).  For  simplicity,  the  depression  storage 
is  omitted  in  this  model,  but  implicitly  is  included  in  the  interception 
storage  capacity  described  in  the  previous  section. 

Thus,  the  water  balance  equation  at  the  ground  surface  is 


I  =  I  -  f. 
e  n  1 


(9) 


in  which  I  is  the  rainfall  excess  rate,  I  is  the  net  rainfall  rate,  and 
e  *  n  ’ 

f.  is  the  infiltration  rate, 
l 

The  ponding  time  and  the  infiltration  rate  are  computed  from  the 
infiltration  model  developed  by  Smith  and  Parlange  (1978).  The  model  was 
developed  from  a  simplified  solution  of  the  equation  for  one  dimensional 
diffusion  of  water  under  gravity,  by  imposing  the  rainfall  pattern  as 
flux-type  boundary  condition.  The  equation  and  the  general  initial  and 
boundary  conditions  are 


90  _  9_ 
at  “  9z 


<D§i) 


(10) 


and, 


t  =  0,  z  >  0,  0  =  0^  ;  (11a) 

0  <  t  <  tp,  z  =  0,  K-D~  =  In(t)  ,  (lib) 


where  0  is  the  soil  water  content;  z  is  the  vertical  distance  from  the  soil 
surface;  D  is  the  moisture  diffusivity;  K  is  the  hydraulic  conductivity;  tp 
is  the  ponding  time;  and  1 n ( t. )  is  time-varying  net  rainfall  intensity 
pattern.  Eq.  (10)  is  also  generally  written  as 


9z  9_  ,_90.  _  dK 
at  90  1  9z'  d0 


(12) 


by  using  an  elementary  identity  of  differentials.  The  primary  assumption 
used  in  the  derivations  is  that  D  varies  rapidly  with  0,  which  implies  that 
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water  flux  within  the  soil  varies  little  with  relative  position. 
Integrating  Eq.  12  with  the  condition  11b,  and  using  an  expression  of 
conservation  of  mass,  Smith  and  Parlange  (1978)  obtained  the  solution 


t  8 

;  I  dt  =  S  (0-0;) 


I  -K(0) 
n 


where  the  subscript  1  refers  to  conditions  at  the  surface,  and  i  refers  to 
initial  conditions.  When  ponding  occurs  8^  =  8g  (saturated  water  content), 
and  the  left-integral  upper  limit  becomes  t  =  t  . 

By  further  assuming  that  K  varies  slowly  in  the  region  near  8s>  or  if 
K  is  much  larger  than  I  ,  Eq.  13  may  be  approximated  by 


Sv  rdt  =  B(0 . )  (I  -  K  ) 
Q  i  np  s 


where  I  is  the  net  rainfall  rate  at  ponding  time,  and  6(8^)  is  a 
parameter  dependent  only  on  soil  type  and  8^.  B  is  theoretically 
identified  as  a  function  of  sorptivity  and  it  can  be  roughly  estimated  from 
B  =  S2/2,  where  S  is  the  soil  sorptivity.  Eq.  14  is  used  in  estimating  the 
ponding  time,  t  .  Where  In(t)  is  a  continuously  differentiable  expression, 
an  analytical  expression  for  t  results  from  Eq.  14.  Engineering  use  of 
this  expression  is  quite  simple.  As  a  rain  pattern  progresses,  the  value 
of  iv* is  calculated  and  compared  with  the  value  of  the  right  hand  side 
of  Eq.  14,  with  I  =  1^ .  This  expression  is  an  inequality  only  up  to  the 
time  at  which  ponding  occurs. 

For  the  same  conditions  on  K  given  above,  the  expression  for  the  decay 
of  infiltration  capacity  for  t  >  t  was  derived  as  was  the  expression  for 
ponding  time.  A  more  general  expression  for  surface  conditions  related  to 
Eq.  (13),  which  reduced  to  Eq.  14,  is 


o  (8-8. )K 
t  .  n 


df  =  J  (0-0. )KdY  S  I  dt 
4*  1  o  n 


Here,  the  variable  of  integration  is  f,  rather  than  8.  This  expression 
applies  for  all  times.  The  value  of  4*  at  the  surface  is  taken  equal  to 
zero  for  t  >  t  ,  although  a  surface  depth  could  be  taken  into  account  if 
desired,  the  effect  of  which  is,  in  fact,  small. 


•  'V* 


Differentiating  Eq.  15  with  respect  to  t  yields 


-  r  IT  °s  (e"ei)  TT^kT5’  ™  =  ?  (Q-ejKdf 

n  dt  f.  1  Un  KJ  4<. 

x  i 


(16) 


Taking  t  >  t^,  the  surface  flux  is  f(t)  <  In(t),  and  integrating  from  t  to 
t  yields 

°  (0-0^  In(f-K)  y  v 

A l  fnr^  *"  'nf-Ttr1  = 

•F  np  np-K 


(t-t  )  J  (e-e.)Kd'f 

P  y  L 

i 


(17) 


Here  A  depends  slightly  on  the  rainfall  rate  pattern.  From  Eq.  14 


A  =  B  =  (I  -K  )  y  1  dt 
np  s  J  n 
r  o 


(18) 


The  function  of  K  in  the  left  integrand  of  Eq.  17  can  be  replaced  by  its 
value  at  saturation,  Kg,  and  combined  with  Eq.  18  to  eliminate  A,  yielding 


t 

Ks  (t-t  )  =  JP  rdt  {~ 

'  "  o  s 


I  -f  I  -K  (I  -K  )f 

_ Q£_s  on  r _ 5£__§ — 1 } 

f-K  K  l(f-K  )I  u 

s  s  np 


(19) 


Eq.  19  is  solved  numerically  to  compute  the  infiltration  decay,  f(t) 
for  t  >  t  .  In  this  model,  Newton's  method  is  applied  and  according  to 
this  method,  the  infiltration  rate  at  the  ith  time  interval  and  kth 


iteration  step  is 


=  f 


F<fi,k-1> 


i.k  ii,k-i  F'd.  ^)  • 


(20) 


The  function  F(f)  and  its  derivative,  F'(f),  are 


I  -f  I  -K 


I  -K 


F(f)  =  K  (t-t  )  -  J*  I  dt  {JE - f-i  tn  (-f-J  7^)1  =  0,  (21) 

0  s  s  np  s 
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and 


To  start  the  computation,  the  following  initial  values  are  assumed: 


f 

0,0 


f(tp)  =  I 


np 


(23a) 


f . 
1,0 


(23b) 


Due  to  the  lack  of  inf i ltrometer  data  in  most  of  the  catchments,  the 

parameter  B  in  Eq.  14  is  presently  estimated  from  the  results  presented  by 

Smith  and  Parlange  (1978),  Fig.  3.  For  simplicity,  the  results  for  Colby 

S.  L.  (swelling)  were  used  in  the  tests  discussed  in  this  report.  With  the 

known  information  of  initial  soil  moisture  deficit,  (0  -0.)  and  the 

si 

saturated  hydraulic  conductivity  Kg  (cm/sec),  the  value  of  B  =  S2/2 

(cm2/sec)  is  estimated  from  the  quasi-linear  functions  plotted  in  Fig.  2. 


2.3  WATER  ROUTING 

The  motion  of  water  moving  either  as  overland  flow  or  channel  flow,  is 
governed  by  the  equations  of  mass  continuity,  momentum  balance,  and  flow 
resistance.  The  flow  routing  scheme  described  in  this  report  is  based  on 
the  kinematic-wave  approximation  to  the  general  equations  of  motion.  The 
following  sections  present  the  basic  assumptions  underlying  the  kinematic- 
wave  theory,  along  with  the  essential  points  of  the  analytical  solution. 
This  solution  is  then  used  to  investigate  the  conditions  under  which 
kinematic  shock  waves  may  be  expected.  Subsequently,  procedures  for  proper 
shock  fitting  are  discussed.  Simplifying  approximations  are  made  which 
allow  closed  form  solutions  and  preserve  the  effects  of  the  shocks.  The 
resulting  approximate  shock  fitting  scheme  is  compared  with  an  existing 
implicit  finite  difference  solution.  The  accuracy  and  efficiency  of  the 
new  scheme  are  illustrated  in  Part  3  by  computing  a  variety  of  unsteady 
flows,  ranging  from  simple  cascades  to  complex  natural  catchments. 
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9  -  9.  INITIAL  SOIL  MOISTURE  DEFICIT 
s  ® 


Fig.  3  Quasi-linear  relation  of  the  second  power  of  soil  sorptivity  to 
initial  moisture  deficit  (after  Smith  and  Parlange,  1978) 


2.3.1  Kinematic  Wave  Approximation 

This  approximation  is  based  on  a  simplification  of  the  Saint-Venant  or 
shallow-water  equations  governing  unsteady  free-surface  flows.  The 
underlying  assumption  is  that  energy  gradients  due  to  local  and  convective 
accelerations  are  negligible  in  comparison  with  gravitational  and 
frictional  effects.  The  momentum  equation  thus  become: 


S0^ 


(24) 


where  and  S^.  are  the  bed  slope  and  friction  slope,  respectively.  The 
continuity  equation  for  water  is  as  usual: 


9a  +  ag  = 

3t  9x 


(25) 


in  which  A  is  the  flow  cross  sectional  area,  Q  is  the  flow  rate  of 
discharge,  x  is  the  downslope  position,  t  is  time,  and  q^  is  the  rate  of 
lateral  inflow  or  outflow  per  unit  length  of  stream. 

Any  suitable  law  of  flow  resistance  can  be  used  to  express  Eq.  24  as  a 
parametric  function  of  the  stream  hydraulic  parameters.  A  widely  used 
expression  is  as  follows: 


Q  =  aAn 


(26) 


where  or  and  n  are  parameters  related  to  channel  (or  plane)  roughness  and 
geometry.  Obviously,  a  and  n  are  functions  of  time  and  position.  The 
space  dependence  can  be  removed  by  making  the  stream  hydraulic  properties 
and  the  lateral  inflow  piece-wise  uniform  in  space  (Fig.  4).  The  flow 
region  is  segmented  into  a  network  of  different  elements  with  properties 
remaining  constant  within  each  element,  but  varying  from  element  to 
element.  This  approach,  known  as  a  kinematic  cascade,  was  introduced  by 
Brakensiek  (1967b)  and  latter  elaborated  upon  by  Kibler  and  Woolhiser 
(1970).  Similarly,  the  time  dependence  of  those  coefficients  and  the 
lateral  discharge  can  be  eliminated  by  assuming  each  piece-wise  constant  in 
time  (Fig.  S).  Then  Eq.  23  and  Eq.  2b  can  be  solved  for  each  cascade 
element  subject  to  given  initial  and  upstream  boundary  conditions  (i.e., 
upstream  inflow  rate,  Qy) -  The  subscripts  U  and  D  will  be  used  to  indicate 
upstream  and  downstream  conditions,  respect i ve 1 y . 
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2.3.2  Analytical  Solution 

Eqns .  25  and  26  can  be  solved  using  the  method  of  characteristics 

(Eagleson,  1970).  The  basic  steps  of  this  method  are  given  below  for 
completeness.  Assuming  that  a  and  n  are  constant  on  a  given  cascade 
element  k  and  combining  Eqns.  25  and  26  yields 


9A 

at 


+  an  A 


n-1  3A 
9x 


(27) 


The  total  differential  of  A  yields 


8A 

at 


dt  +  dx  = 
dx 


dA 


(28) 


Since  the  derivatives  in  Eqns.  25,  27,  28  do  not  exist  along  the 
characteristic  paths,  the  determinant  of  the  coefficient  matrix 
corresponding  to  those  equations  must  vanish  at  points  on  the 
characteristics.  This  requirement  thus  yields 


dx 

dt 


=  omA 


n-1 


1  (1  -  -) 
na11  Q  n 


(29a) 


or, 


dt 

dx 


anA 


n-1 


Obtaining  the  invariants  of  the  solution  gives 


dA  _ 
dt  q£ 


(29b) 


(30a) 


or, 


dA 

dx 


anA 


n-1 


(30b) 


The  sets  of  Eqns.  29A,  30A  or  29b,  30b  represent  the  characteristic 
form  of  the  solution  to  the  kinematic-wave  approximation.  They  show  that 
kinematic  waves  possess  only  one  system  of  characteristics .  Accordingly, 
kinematic-wave  routing  cannot  be  used  in  situations  where  there  are 
downstream  flow  controls.  Integrating  Eqns.  30a  and  30b  with  the  initial 
and  boundary  conditions 
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Aj  =  A(x,ty) ,  tQ  =  0, 

(31) 

Ay  —  A(Xy,t),  xQ  -  0, 

(32) 

results  in  the  following  respective  expressions  for  the  flow  conditions 
along  any  characteristic  t  =  t(x)  (Harley  et  al.,  1970): 


and 


A(x,t)  =  A  +  f  q. { x ,  £(x) ]  d£  , 
l0 


A(x,t)  =  JAo  +  h  f  do)  n  , 


(33a) 


(33b) 


or 


Q(x»t)  =  aA[j  +  J  q£(n,t(n)l  dr). 

xo 


(33c) 


Substituting  Eqns.  33a  and  33b  into  Eqns.  29a  and  29b,  respectively,  and  on 
integration  yield 


t  t' 

x-xQ  =  an  /  {/  q£[x,£(x)J  d£  +  AQ} 

C0  t0 


n- 1 


dt', 


(34a) 


and 


1-n 


t_t0  =  k  f  {a  f  q£ln’t(n)J  dn  +  aJ}  n  dx’  . 
xo  x0 

In  these  equations,  AQ  represents  the  initial  flow  area,  Aj 


(34b) 


along  the 

characteristic  (shown  in  Fig.  5),  and  the  upstream  inflow  area,  Ay  = 
1 1/n 


[Qu(t)/o] 


along  characteristics  like  Cj,  C£,  and  C^.  The  above 


integrals  are  functionally  integrable  and  yield  simple  expressions  when 
q£(x,t)  is  either  an  explicit  function  of  x  and  independent  of  time,  or 
vice  versa.  In  most  applications,  only  discrete  distributions  of  lateral 
inflows  are  available.  For  instance,  rainfall  intensities  and/or 
infiltration  rates  are  usually  regarded  as  lateral  inflows  in  hydrologic 
models.  In  these  cases,  it  has  been  a  common  practice  to  use  lateral 
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inflows  having  a  piece-wise  uniform  distribution  in  space  (Fig.  4)  and  a 
piece-wise  constant  variation  in  time  (Fig.  5).  This  form  of  lateral 
inflow  distribution  is  also  adopted  in  this  paper.  Thus,  assuming  q^(x,t) 
remains  constant  within  small  space  and  time  increments  the  above 
expressions  may  be  explicitly  integrated.  Discretizing  Eqns .  33  and  34 
between  two  points  and  (x^tj)  on  a  characteristic  path  (Fig. 

5)  yields 


A.  .  =  A .  _  .  ,  +  qn  .  At .  , 

i.J  i-i»j-i  £,J  J 

1 

A.  .  =  (A.n.  .  +  Ax .  )n  , 

i.J  i-l, J-l  a  i 


Qi,j  =  Qi-l,j-l  +  **i  ’ 


and 


Ax .  = 
i 


[  (q„  .  At .  +  A .  .  )n  -  A . n  .  ,  J  , 

q„  .  £,J  J  x-1, j-l  i-l,  j-l 

* » J 


At.  =  -1—  [(^J-  Ax.  +  A,n  .  _)"  -  A.  ,  .  J  , 
j  a  i  x-i, j-r  ’ 


1 

,n 


(33'a) 
(33'b) 
(33'c) 
(34 'a) 

(34 'b) 


where  Ax  =  x.  -  x.  ,,  At.  =  t.  -  t.  ,  and  q»  .  is  the  rate  of  lateral 

inflow  assumed  constant  over  the  time  step  At.  and  the  space  increment  Ax.. 

J  1 
The  last  two  equations  are  not  defined  when  q„  .  =0;  in  this  case  Eqns. 

“ ,  J 

29a  and  29b  yield 


a  »n-l 

Ax.  =  an  A.  ,  .  ,  At . 
i  i-l, j-l  j 


and 


A.  1  ,1-n  . 

At .  =  —  A .  .  .  .  Ax . 
j  an  i-l, j-l  i 


134’ c) 


(34' d) 


Eqns.  34'  are  used  to  trace  the  characteristic  path  by  considering 

either  Ax^  or  At^  as  a  dependent  variable  and  choosing  a  suitable  value  for 

the  other.  This  increment  must  be  chosen  so  that  the  lateral  inflow  can  be 

assumed  constant  within  the  intervals  At.  and  Ax..  Eqns.  33'  are  used  in 

J  i 

turn  to  compute  the  flow  conditions  on  the  same  characteristic.  Eqns.  33' 
and  34'  are  explicit  and  independent  of  each  other,  and  therefore,  they  may 
be  used  in  several  different  combinations.  In  this  paper,  the  following 
scheme  is  used  for  computational  efficiency.  Since  the  lateral  inflow  is 
uniform  along  a  cascade  element  and  piece-wise  constant  with  time,  the  best 
choice  is  to  take  the  time  increment  At.  as  the  independent  variable.  This 
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time- increment  is  used  in  Eq.  33'a  to  compute  A^  ^ ,  and  Eq.  26  is  used  to 

obtain  Q.  ..  Then  Eq.  33'c  is  solved  for  Ax..  These  steps  constitute  the 
1  >  J  i 

scheme  used  to  trace  a  characteristic  across  the  cascade  element  and  to 

compute  the  flow  conditions  along  its  path. 

In  general,  a  characteristic  does  not  intersect  the  downstream 
boundary  exactly  at  the  beginning  or  end  of  a  time  step.  A  modification  of 
the  proposed  scheme  is  used  to  determine  the  time  of  intersection.  The 
last  space  increment  Ax^  =  x^_j,  where  x^  is  the  length  of  the  cascade 
element,  is  substituted  in  Eq.  33'c  to  compute  and  then  Eq.  26  is  solved 
for  Ap.  This  value  is  then  introduced  into  Eq.  33'a  to  solve  for  the  time 
increment,  At..,  and  the  time  of  intersection  t_  =  t„  ,  +  AtVI.  The 

discharges  existing  on  all  the  characteristic  paths  at  the  time  of  their 
arrival  to  the  downstream  boundary  define  a  discrete  outflow  distribution. 
This  discrete  distribution  is  smooted  out  by  interpolation  to  obtain  a 
continuous  outflow  hydrograph.  The  problem  of  hydrograph  computation  for 
any  given  cascade  element  is,  thus,  completely  specified  once  the  initial 
flow  areas  along  the  element  at  time  zero  (Eq.  31)  and  the  inflow 
hydrograph,  coming  from  the  upstream  element  (Eq.  32),  are  known. 


2.3.3  Formation  of  Shock  Waves 

The  foregoing  analytical  solution  remains  valid  as  long  as  the 
characteristic  paths  do  not  intersect  each  other.  When  this  occurs,  the 
preceding  theory  does  not  give  a  unique  solution  since  there  is  more  than 
one  characteristic  passing  through  the  intersection  point,  causing  flow 
discontinuity.  These  discontinuities  have  properties  analogous  to  those  of 
shocks  that  arise  in  gas  dynamics  theory  and,  by  analogy,  are  generally 
known  as  kinematic  shock  waves.  Eq.  29a  shows  that  the  celerity  of  a  wave 
moving  along  a  characteristic  path  is  a  function  of  the  flow  area  A.  This 
dependence  on  A  produces  a  nonlinear  distortion  of  the  wave  as  it 
propagates.  Waves  with  higher  values  of  A  travel  faster  and  finally 
overtake  lower  ones  leading  to  the  breaking  of  the  wave  profile  as 
illustrated  in  Fig.  6.  In  this  figure  the  lateral  inflow  has  been  ignored 
for  simplicity  and,  thus,  the  flow  areas  do  not  change  along  the 
characteristic  paths  (Eq.  30).  The  wave  profile  at  any  instant  t  =  tj  >  0 
is  obtained  by  projecting  the  corresponding  values  of  A  on  the  verticals 
passing  through  the  intersections  of  the  characteristic  paths  and  the  line 
t  =  t,.  The  wave  breaking  begins  at  the  time  t  =  t_  when  the  free  surface 
profile  first  develops  an  infinite  slope.  In  the  x-t  plane  this  breaking 
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coincides  with  the  intersection  of  two  characteristics,  where  a  single¬ 
valued  flow  area  no  longer  exists.  This  sudden  change  in  flow  area  is 
interpreted  as  the  initiation  of  a  shock.  After  its  formation,  the  shock 
proceeds  downstream  with  its  own  velocity,  as  discussed  in  the  next 
section. 

However,  a  shock  will  be  formed  only  if,  given  any  two  consecutive 
characteristics,  the  first  one  is  less  steep  than  the  following  one.  For 
example,  Fig.  5  shows  the  characteristics  and  passing  through  the 
points  E  and  D  at  the  same  time  t'.  Then,  the  necessary  condition  for 
shock  formation  is  that 


dx  (1) 
f— 1 
dt  E 


[  — ] 
ldtJ 


(2) 

D 


(31) 


where  the  superscripts  refer  to  the  consecutive  characteristics  and  Ca¬ 
using  Eqns.  26,  29a  and  33a,  Eq.  35  can  be  expressed  as 


1 


1 


1 


r_  (2) ,nk  f.  (1) ,nk  .  ,  "k  A. 

lQU,k  1  ‘  lQU,k  1  (V  q2,k  Atj  +  1 


(36) 


in  which  the  subscript  k  indicates  flow  variables  associated  with  the  kth 

cascade  element  and  At.^,  =  t.^,  -  t..  This  inequality  indicates  that 

J+1  J+1  J 

characteristics  emanating  from  the  line  t  =  0  cannot  intersect  for  either 
dry  or  uniform  initial  conditions.  It  is  also  clear  from  Eq.  36  that 
characteristics  originating  during  steady  upstream  inflow  or  along  the 
falling  limb  of  the  upstream  hydrograph  (i.e.,  characteristics  in  Fig. 
5)  will  not  intersect,  unless  the  flow  law  changes  from  turbulent  to 
laminar.  These  two  zones  are  then  free  from  shock  formation.  This  does 
not  mean  that  shocks  originating  upstream  may  not  propagate  into  these 
zones.  On  the  other  hand,  the  characteristics  originating  along  the  rising 
limb  of  the  upstream  hydrograph  may  intersect  depending  upon  the  relative 
magnitude  of  the  variables  appearing  in  Eq.  36.  For  this  reason,  the 
domain  bounded  by  the  characteristics  emanating  from  boMi  ends  of  the 
rising  limb  can  be  regarded  as  the  probable  shock-forming  zone.  The 
solution  in  the  shock-forming  zone  is  influenced  by  the  upstream  boundary 
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conditions  of  the  kth-element  which,  in  turn,  depend  on  the  initial 
conditions  on  the  element  k-1.  For  example,  applying  Eq.  36  to  the 
characteristics  passing  through  the  points  M  and  N  indicated  in  Fig.  5,  and 
using  Eqns.  29a  and  33c,  one  obtains 


k-1 

(JLJ.) 

v  /M  ' 


1. 

n, 


k-1 


’t.k'Vi 


(37) 


This  inequality  shows  that  for  uniform  topography  (a,  =  a  ,  n,  =  n  )  a 

K“ 1  K  K" i  K 

shock  will  be  formed  in  element  k  whenever  q^  k  ^  >q^  whereas,  under 
conditions  of  uniformly  distributed  lateral  inflow  (q0  =  q  ),  a  shock 

will  occur  wherever  akl  >  ak  for  nkl  =  nk  and  nk  J  >  nk  for  ak  =  ak* 

Harley  et  al.  (1970).,  Kibler  and  Woolhiser  (1970),  and  Li  et  al. 

(1976a, b)  also  studied  the  foregoing  necessary  conditions.  The  latter  also 
discussed  the  sufficient  conditions  for  two  given  characteristics  to 
intersect  within  the  boundaries  of  the  cascade  element.  They  suggested 
taking  two  consecutive  characteristics  far  apart  enough  so  as  to  avoid 
their  intersection  within  the  element  (i.e.,  characteristics  C'  and  C"  in 
Fig.  5).  However,  once  the  necessary  condition  is  met,  the  shock  wave  S 

may  form  within  the  cascade  element  and  travel  all  the  way  to  the 

downstream  end  (this  is  discussed  in  the  next  section).  Therefore,  not 
considering  the  solution  domain  between  C'  and  C"  serves  only  to  ignore  the 
possible  existence  of  the  shock  but  does  not  really  avoid  it.  Given  that 
the  shock  is  an  intrinsic  part  of  the  solution  to  the  kinematic-wave 
approximation,  a  better  routing  procedure  is  one  that  considers  the 
presence  of  the  shock  and  its  effect  cn  the  outflow  hydrograph.  Such  an 
approach  is  discussed  in  the  next  section. 


2.3.4  Approximate  Solution  With  Shock  Fitting 

After  breaking  occurs  the  kinematic  wave  theory  ceases  to  be  a 
strictly  valid  description  of  the  physical  process,  because  the  flow  area 
is  inherently  single-valued.  However,  the  foregoing  formulation  can  still 
be  utilized  by  allowing  discontinuities  in  the  solution.  Lighthill  and 
Whitham  (1955)  proposed  replacing  the  shock,  as  a  first  approximation,  by  a 
discontinuous  wave  (BC,  Fig.  6)  that  produces  the  appropriate  variations  in 
discharge  and  flow  area  as  it  moves  downstream.  These  authors  derived  the 
velocity  of  the  shock  by  considering  the  continuity  of  flow  and  the  rate  of 
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discharge  crossing  the  shock  front.  Applying  their  result  to  an  arbitrary 
points  (x,t)  on  the  shock  path  resulted  in 


u(x  t)  =  dx  Q  Jx2t) 

*•  ’  J  dt  ,b 


Q  (x,t) 


(38) 


A  (x,t)  -  A  (x , t) 


where  the  superscripts  b  and  a  indicate  flow  conditions  behind  and  ahead  of 
the  shock,  respectively.  While  the  shock  is  moving  downstream, 
infinitesimal  waves  traveling  along  characteristic  paths  will  join  the 
shock  from  ahead  and  behind  continuously  modifying  the  shock  strength  and 
velocity.  Therefore,  the  locus  of  all  the  points,  where  those 
characteristics  intersect  each  other,  defines  the  path  of  the  associates 
shock. 

The  important  steps  involved  in  routing  a  shock  include  finding  the 
position  where  it  originates  and  then  tracing  the  shock  downstream.  An 
analytical  method  for  finding  the  position  where  the  shock  originates  was 
discussed  by  Whitham  (1974).  Because  of  its  computational  intricateness, 
his  approach  was  not  used  in  the  present  study,  but,  instead,  the  following 
simpler  technique  was  developed.  It  is  clear  from  Eq.  36  that  once  it  is 
satisfied,  the  characteristics  tend  to  converge.  Thus,  a  shock  may 
originate  within  the  same  element  or  in  one  of  the  following  cascade 
elements,  depending  upon  the  magnitudes  of  a^,  n^,  and  q^  ^  as  they  appear 
in  Eq.  .36.  In  this  paper,  the  problem  of  finding  the  exact  location  of  the 
shock  origin  is  avoided  by  discretizing  the  upstream  inflow  hydrograph.  By 
doing  this,  small  artificial  shocks  are  introduced  all  along  the  rising 
limb  of  the  hydrograph,  but  they  are  routed  across  the  characteristic  plane 
only  in  those  zones  where  Eq.  36  is  satisfied.  The  remainder  of  the  inflow 
hydrograph  is  routed  according  to  the  continuous  solution  from  the 
characteristics.  Thus,  the  shock  fitting  approximations  introduced  in  this 
paper  influence  only  those  zones  where  shocks  form,  whereas  the  accuracy  of 
the  analytical  solution  is  preserved  in  the  other  zones.  The  procedure  is 
illustrated  in  Fig.  7,  where  individual  characteristics  are  traced  one  at  a 
time,  starting  at  the  upstream  boundary  of  the  characteristic  plane.  They 
are  assumed  to  emanate  from  the  half-time  levels  t^,  tj+^,  t2+V  etc'* 
where  the  continuous  and  discretized  upstream-hydrographs  coincide.  The 
shock-forming  condition,  Eq.  36,  is  tested  at  each  new  time  step.  The 
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discharges  and  Qy^  appearing  in  this  relationship  are  replaced  by  the 

discretized  inflow  rates  at  the  new  and  previous  time  steps,  respectively. 

In  addition,  the  product  ^At^  +  ^  is  replaced  by  the  integral  of  the 

lateral  inflow  hydrograph  over  the  interval  t .  ,  i  t  $  t .  ,  ,  to  account  for 

J-"!  J  +  "5 

the  fact  that  the  characteristic  paths  arise  from  t.  ,  instead  of  t. 

J+!s  J  +  l 

Hence,  if  j  =  2,  for  instance,  Eq.  .36  is  given  by 


i_  J_  I_ 

(QU^)0"  -  (Qy,)^  >  («/“  F2, 

where  F2>  Qy2>  an<*  Qy3  are  defined  in  Fig.  7.  When  Eq.  36  is  not 
satisfied,  the  corresponding  characteristics  are  routed  one  at  a  time  using 
Eqns.  33  and  34,  as  illustrated  by  C^,  C ^ ,  and  C  in  Fig.  7.  On  the  other 
hand,  if  Eq.  36  is  met,  an  artificial  shock  is  introduced  at  the  time  level 
t ^  and  then  routed  using  the  scheme  presented  below. 

The  problem  of  tracing  the  shocks  in  the  characteristic  plane  has  been 
discussed  in  detail  by  Whitham  (1974).  It  consists  essentially  of 
determining  the  characteristics  intersecting  each  other  on  the  shock  path, 
and  the  location  of  the  point  of  intersection.  The  solution  to  this  so 
called  three-point  problem,  requires  the  simultaneous  solution  of  the 
characteristic  Eqns.  34  and  the  expression  for  the  shock  velocity,  Eq .  38. 
This  method  was  used  by  Kibler  and  Woolhiser  (1970),  who  developed  an 
iterative  scheme  to  route  shocks  generated  on  a  three-plane  cascade  under 
constant  lateral  inflow.  In  principle,  this  approach  can  be  extended  to 
include  more  complex  geometries  with  varying  lateral  inflow.  However,  the 
computational  details  become  excessively  complicated,  particularly  when 
shock  interactions  occur.  Therefore,  an  alternate  approach  is  developed  in 
this  paper  by  restricting  the  above  artificial  shocks  to  small  amplitudes 
(weak  shocks).  Consider  the  characteristic  pair  C^,  intersecting  at  the 
point  P  on  the  shock  path  (Fig.  7),  and  a  similar  pair  C^,  joining  the 
shock  at  an  arbitrary  point  B.  If  the  shock  is  weak,  the  flow  areas 
carried  by  and  (ahead  of  the  shock)  will  not  differ  significantly 
from  each  ether.  Similarly,  and  will  have  approximately  equal  flow 
areas  behind  the  shock.  Since  B  is  arbitrary,  and  Clj  may  represent  any 
of  the  infinite  pairs  of  characteristics  joining  the  shock  (behind  and 
ahead  of  it)  along  its  trajectory.  It  is  thus  reasonable  to  postulate  that 
the  intersecting  characteristics  defining  the  shock  path  do  not  deviate 
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significantly  from  it.  Then,  as  a  further  expedience,  we  assume  that  the 
shock  path  is  defined  by  two  sufficiently  close  (overlapping) 
characteristics,  one  representing  the  characteristics  ahead  of  the  shock 
and  the  other  representing  the  characteristics  behind  the  shock.  These  two 
characteristics  will  carry  flow  areas  (depths)  ahead  and  behind  the  shock 
as  given  by  Eq.  33a  or  33b.  The  common  path  of  these  characteristics,  as 
well  as  the  shock,  is  obtained  by  substituting  Eqns .  33a  and  26  or  Eqns. 
33b  and  33c  into  Eq.  38;  after  integration  they  yield 

t  ,  t ' 

X-X0  =  ~ET~a  $  *lA0  +  q£(x»Ux))  d£)n  "  IAq  + 

VA0  t0  t0 


/  q„(x>4(x))  d£]n}  dt'  , 


(39a) 


XX  — 

t-tn  =  TT-^T  /  { I (A«)n  +  l  S  q0(n,t(n))  dp]"  -  [(A*)“  + 


xo 


l  •  /  q£(n,t(n))  dnln]  dx’  , 
xo 


(39b) 


where  (A^,  Qq)  and  (A^,  Q^)  are  the  known  flow  conditions  ahead  and  behind 
the  shock,  respectively,  at  a  given  point  (Xg,t^)  on  the  shock  trajectory. 

Consider  now  any  two  consecutive  points  B(x^  ^.t^^)  and  D(x^,t^)  on  a 
shock  trajectory  S  (Fig.  7).  Applying  Eqns.  39  between  these  two  points 
and  discretizing  the  results,  gives 


Ax.  = 


(n+l)q. 


uvtrAh.i-il"n 


'V'-j  + 


At.  = 


J  (n+l)q£  (Q 


an  i,  i  ,  ,  .,b  ,i 

b - - :  u<r  *  (Ai-,,j-,) 

i-1 ,  j-1 


a  >  /Aa  ,n,  n  ,Ab  .n+1 

lT  **i  +  (Ai-lfj-i)  1  -  (Ai-i,j-i)  + 


(Aa  )n+1} 


(40b) 


in  which  the  lateral  inflow  is  taken  as  constant  over  the  finite  increments 


Ax.  and  At..  Eqns .  40  are  not  valid  when  q.  =  0.  In  this  case,  a  discrete 
i  J  * 

approximation  of  Eq.  38  is  used.  Similar  to  the  characteristic  wave 
routing,  At^  is  taken  as  the  independent  variable.  Eq.  33 'a  is  used  to 

compute  A3  .  and  Ab  .  and  using  these  two  in  Eq.  26  0a  .  and  Q b  .  are 

^*.1  i  i  J  i )  J 

determined.  Then  Eq.  40a  is  used  to  compute  Ax^.  These  steps  constitute 

an  explicit  scheme  that  is  used  to  route  the  shock  like  a  characteristic 

wave.  Similar  to  the  characteristic  wave,  the  scheme  is  modified  near  the 

downstream  boundary.  Here,  Eqns.  33' c  and  26  are  used  to  compute  the  flow 

conditions  ahead  and  behind  the  shock,  and  Eq .  40b  is  used  to  find  the 

exact  arrival  time  of  the  shock.  Although  this  scheme  is  based  on  an 

approximation  that  deviates  from  the  uniqueness  of  the  solution,  its  use 

does  not  significantly  affect  the  accuracy  of  the  computations.  In  the 

foregoing  discussion  the  shock  fronts  have  been  considered  mathematically 

as  flow  discontinuities  fitted  into  regions  where  multiple  values  of  the 

solution  exist.  Physically  they  will  not  be  discontinuities ;  instead,  the 

fronts  will  have  finite  lengths  induced  by  diffusive  effects  as  well  as  by 

breaking  of  the  waves.  In  many  instances,  the  front  will  take  the  smooth 

shape  of  a  monoclinal  flood  wave  (Whitham,  1974).  Because  of  these 

smoothing  tendencies,  and  for  the  sake  of  computational  expedience,  the 

shock  front  can  be  approximated  by  a  linear  profile  as  shown  in  Fig.  7 

(segment  EF).  Thus,  the  flow  condition  at  the  shock  itself  can  be  computed 

as 

A(x , t)  =  |  (Aa(x,t)  +  Ab(x,t)].  (41) 

This  could  be  viewed  as  a  method  to  make  an  approximately  continuous 
solution  from  one  which  was  made  discontinuous  by  discretizing  the  inflow 
hydrograph.  The  preceding  routing  approach,  called  herein  a  propagating 
shock-fitting  (PSF)  scheme,  permits  the  use  of  essentially  the  same 
numerical  procedure  to  route  both  characteristic  and  shock  waves.  This 
scheme  is  particularly  efficient  when  the  initial  flow  conditions  are 
either  dry  or  uniform,  and  only  the  outflow  hydrograph  is  of  interest.  In 
these  cases  the  calculations  are  performed  only  in  the  x-t  subdomain 
bounded  by  the  lines  C.  and  C  (Fig.  7)  where  C  is  the  characteristic 
reaching  the  downstream  boundary  at  the  end  of  the  routing  interval.  All 
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the  characteristics  emanating  from  the  line  t  =  0  are  parallel  to  C^. 
Hence,  the  outflow  conditions  at  times  t^,  t^,  t^,...  are  equal  to  those 
computed  at  the  points  1,  2,  3,...  on  C^.  Within  the  above  subdomain  the 
step-by-step  integrations  along  the  wave  trajectories  use  simple  algebraic 
relationships.  Moreover,  values  computed  at  interior  points  do  not  have  to 
be  stored  since  only  outflow  conditions  are  needed. 

In  the  same  manner  that  shocks  arise  from  the  intersection  of 
characteristic  waves,  they  can  also  meet  with  other  shocks  to  form  new 
shocks.  In  addition,  shocks  introduced  in  the  shock-forming  zone  of  a 
given  cascade  element  will  propagate  into  the  downstream  elements  of  the 
cascade  interacting  with  each  other  and  creating  new  shocks.  These  shock 
interactions  cannot  be  treated  by  the  above  PSF  scheme,  because  it  tracks 
only  one  wave  at  a  time.  For  this  reason,  a  further  simplification  is 
introduced  which  consists  of  restricting  the  shock  interactions  to  the 
junctions  of  the  kinematic  cascade.  The  shocks  emanating  from  the 
discretized  inflow  hydrograph  are  tracked  across  the  shock  forming  zone 
using  the  explicit  scheme,  Eqns .  40.  When  all  the  shocks  have  been 
projected  to  the  downstream  boundary,  their  fronts  are  smoothed  out  using 
Eq.  41.  A  smooth  outflow  hydrograph,  incorporating  the  overall  effect  of 
the  shocks  formed  upstream  is  thus  obtained.  The  interaction  of  shocks 
carried  by  converging  outflow  hydrographs  is  then  simulated  by  simple 
superposition  of  these  hydrographs.  The  resulting  outflow  hydrograph, 
which  satisfies  flow  continuity,  is  in  turn  used  as  upstream  boundary 
condition  to  the  next  element  and  the  procedure  is  repeated.  This  method 
will  be  called  an  approximate  shock  fitting  (ASF)  scheme.  This  scheme 
computes  the  outflow  hydrograph  at  any  junction  of  the  cascade  reflecting 
the  overall  effects  of  the  shocks  formed  in  the  upstream  elements. 

2.4  SEDIMENT  ROUTING 

Sediment  yield  from  agricultural  catchments  is  controlled  by  physical 
principles  governing  the  detachment  and  transport  of  sediment  particles. 
The  source  of  sediment  erosion,  or  sediment  supply,  is  the  detachment  by 
raindrop  impact  and  by  runoff.  If  the  sediment  load  from  upstream  areas  is 
less  than  the  potential  transport  capacity  of  the  flow,  the  supply  is 
depleted  and  erosion  occurs.  If  the  sediment  load  is  greater,  deposition 
occurs.  These  processes  are  all  interrelated  and  must  satisfy,  locally, 
the  conservation  principle  of  sediment  mass.  Therefore,  in  addition  to  the 
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equations  of  continuity  and  momentum  for  water,  additional  equations  are 
required  for  sediment  routing.  These  are  the  sediment  continuity  equation, 
and  relations  describing  sediment  supply  and  transport.  These  equations 
are  presented  below,  and  they  are  equally  applied  to  overland  and  channel 
flow. 


2.4.1  Sediment  Continuity  Equation 

The  continuity  equation  for  the  size  classes  k  =  1,  2,  ...,  n,  forming 
the  sediment  load  can  be  written  as 
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where  G  .  is  the  sediment  transport  rate  by  volume  per  unit  time,  X  is  the 

S  K 

soil  porosity,  P  is  the  wetted  perimeter,  z^  is  the  land  surface  elevation, 
g^k  is  the  lateral  inflow,  and 

Ck  =  Gsk/Q  <M> 

is  the  fraction  concentration  by  volume  (Fig.  8a).  The  concentration  and 
bed  elevation  for  the  total  load  are 
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the  third  term  in  Eq.  42  will  be  denoted  by 
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This  term  can  be  envisioned  as  a  sediment  supply  term  for  exchange  between 
the  flow  and  the  detached  bed  material  during  erosion  or  deposition. 
Assuming  that  the  sediment  moves  essentially  at  the  same  average  velocity 
of  flow,  V,  Eq.  42  can  be  expressed  as 
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Taking  V  as  approximately  constant  over  small  space  and  time  intervals, 
yields 


3A  .  3A  . 
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(48) 


in  which  A  ,  =  C,  A  =  G  . /V  is  the  volume  of  sediment  fraction  present  in 

sk  k  sk 

the  flow  per  unit  length.  Eq.  48  is  used  to  track  the  aggradation  or 
degradation  of  the  bed  as  explained  in  section  2.4.4. 

The  present  model  uses  a  power  function  to  relate  the  wetted  perimeter 
to  the  flow  area.  That  is 

P  =  aAb 

where  a  and  b  are  coefficients  depending  on  the  cross-sectional  shape  of 
the  flow  reach.  In  particular  a  =  1.0  and  b  =  0.0  for  overland  segments. 

2.4.2  Sediment  Transport  Formulas 

These  formulas  are  used  to  determine  the  sediment  transport  capacity 
for  a  specific  set  of  flow  and  sediment  characteristics.  The  formulas  used 
in  the  present  model  were  selected  after  a  recent  study  by  Alonso  et  al. 
(1981).  They  compared  the  predictions  of  nine  transport  formulas  with 
flume  and  field  data.  The  comparison  was  based  on  40  field  measurements, 
523  flume  experiments,  and  176  tests  on  concave  slopes,  with  sediments 
ranging  from  coarse  sands  to  very  fine  soil  particles.  As  expected,  no 
formula  satisfactorily  represented  the  entire  spectrum  of  sediment  and  flow 
characteristics.  Nevertheless,  three  of  the  tested  formulas  gave 

satisfactory  estimates  of  transport  capacity  over  different  subsets  of  the 
data  range. 

The  total  load  formula  of  Yang  (1973)  best  estimated  streamflow 

carrying  capacity  in  the  range  of  fine  to  coarse  sands.  The  total  load 
formula  of  Laursen  (1958)  predicted  reasonably  well  small  streamflows 
carrying  very  fine  sands  and  silts.  It  should  be  used  with  some 

reservations,  however,  for  computing  transport  of  lighter  materials  such  as 
soil  aggregates.  The  bed  load  formula  of  Yalin  (1963)  can  be  used  to 

compute  sediment  transport  capacities  for  overland  flows.  These 
conclusions  are  graphically  summarized  in  Fig.  9.  Each  of  these  formulas 
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is  presented  below  as  the  originator  intended,  but  rewritten  in  terms  of 
dimensionless  parameters,  (Alonso  et  al.,  1981).  Where  the  formulations 
required  graphical  solutions  (i.e.,  determination  of  threshold  conditions 
from  Shield's  curve),  analytical  equivalents  (not  shown  here)  have  been 
worked  out  to  facilitate  their  use  in  digital  computation. 

(i)  Total  load  formula  of  Yang  (1973):  Yang  based  his  formula  on  the 
premise  that  total  load  is  dominated  by  the  rate  of  potential  energy 
expenditure  per  unit  weight  of  water.  He  used  this  concept  and  dimensional 
analysis  to  derive  his  formula,  the  coefficients  of  which  were  determined 
from  a  large  set  of  laboratory  data.  This  formula  is: 

“’k  =  6k  Zk  (V/U*)(10<1>"6/S]’  (49) 

where 

«J»  =  5.435  -  0.286  log  (wkdk/v)-0 . 457  log  (u*/wk)  + 

(1.799  -  0.409  log  (wkdk/v)  -  0.314  log  (u*/wk) ] 

log  (VSQ/wk  -  VcSQ/wk),  (50) 

Vc/w  =  2.5/[log(u*dk/v)  -  0.06]  +  0.06,  0<u*dR/v<70,  (51a) 

V  /w  =  2.05,  u.d./v  >  70.  (51b) 

c  w  k 

(ii)  Total  Load  Formula  of  Laursen  (1958):  based  mostly  on  heuristic 
considerations ,  Laursen  developed  a  formula  relating  the  load  concentration 
to  the  relative  roughness  and  excess  tractive  force.  This  relationship  is 
corrected  by  an  empirical  function  of  (uj,./w)  which  accounts  for  the 
effectiveness  of  turbulence  in  suspending  the  bed  material.  The 
contributions  of  each  size  fraction  are  added  up  to  yield  the  total  load. 
The  formula  is: 
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Fig.  9  Range  of  applicability  of  sediment  transport  formulas 


where  f(u^/w^)  is  the  empirical  function  obtained  by  Laurst n  from  t  1  um< 
experiments . 

(iii)  Bed  Load  Function  of  Yalin  (1963):  this  formula  is  based  on  a 
theoretical  analysis  of  saltating  particles,  in  which  it  is  proposed  that 
the  bed  load  rate  is  related  to  the  range  of  the  particles'  saltation 
rather  than  to  the  number  of  particles  participating  in  the  motion  fhr 
resulting  formula  is  of  the  excess-shear  type,  and  was  calibrated  on  a 
limited  set  of  laboratory  data.  This  formula  has  the  form 

*k  =  0.635  0 1'2  II  -  (0ck/ek)l  |(i/eck) 

-  (1/ao  0cR)  ln(  1  +  ao]  ,  (5.3) 


in  which 


a  =  2.45  S"2/S  0^.  ,  (54) 

ck 

a  =  (e5O/0ck)  _  1  =  fu*/eck(s‘i)  gdk]  ' 

In  Eqns.  49  through  55  <t>k  is  the  dimensionless  volume  transport  rate, 
0^  and  are  the  mobility  number  and  relative  roughness  based  on  the  d^ 
fraction  size,  w^  is  the  settling  velocity  of  sediment,  and  the  subscript  c 
denotes  threshold  conditions.  These  parameters  are  defined  as 


\  =  (Gsk/yPSi/Us-Dgdj^, 

(r.b) 

6k  =  f  Cs-i)gdk] , 

(57) 

Zk  =  y/dk> 

(58.) 

where  S  is  the  specific  gravity  of  the  sediment  fraction,  u.v  is  the  bed 
shear  velocity,  and  y  and  V  are  the  specific  weight  and  kinematic  viscosity 
of  water,  respectively. 
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2.4.3  Equations  of  Sediment  Supply 

Soil  eroded  from  land  areas  with  distributed  and  concentrated  flows  is 
the  source  of  most  of  the  sediment  transported  in  a  catchment  system.  Soil 
erosion  is  a  complex  process  of  detaching  soil  particles  and  transportng 
them  downslope  through  the  action  of  raindrop  impact  and  runoff.  Erosion 
begins  when  raindrops  strike  the  land  surface  and  detach  soil  particles  by 
splash.  Whenever  the  soil  surface  is  not  protected  by  vegetation,  or  any 
other  cover,  raindrops  can  detach  very  large  quantities  of  soil.  Most  of 
this  eroded  soil  is  moved  downstream  by  surface  runoff.  When  runoff 
reaches  sufficient  intensity,  the  rate  of  soil  erosion  becomes  also 
dependent  on  the  runoff  characteristics  and  on  the  susceptibility  of  the 
soil  to  the  forces  of  the  flowing  water.  The  following  sections  describe 
the  approaches  used  in  the  present  model  to  simulate  these  major  sources  of 
eroded  sediment. 


2.4.3. 1  Soil  Detachment  by  Raindrop  Impact:  Past  research  on  the  process 
of  soil  detachment  by  raindrop  impact  suggested  that  the  rate  of  detachment 
is  proportional  to  the  square  of  the  net  rainfall  intensity  (Meyer  and 
Wischmeier,  1969).  However,  bare  soils  differ  in  their  susceptibility  to 
detachment  by  raindrop  impact  due  to  a  variety  of  properties  such  as 
primary  particle  size  distribution,  soil  structure,  organic  matter  content, 
etc.  Therefore,  an  erodibility  parameter,  a^,  is  introduced  to  describe 
the  rate  of  erosion  as 

D  =  a  I2.  (59) 
r  r 

This  relationship  has  been  recently  cor.i.rmed  by  Meyer  (1980)  who  reported 
extensive  field  studies  showing  that  Eq.  59  is  a  valid  approximation  for  a 
wide  range  of  soils  and  cropping  conditions.  Eq.  59  is  the  basic  equation 
for  rainfall  detachment.  Parameters  must  be  added  to  account  for  other 
factors  that  influence  the  rate  of  detachment.  Meyer  (1980)  has  shown  that 
the  I2  law  is  not  affected  by  the  stage  of  canopy  development,  but  the  af 
coefficient  is  dependent  on  the  erodibility  of  the  soil  and  decreases  as 
the  vegetative  cover  changes  from  first  growth  to  full  canopy.  To  account 
for  this  effect  a  cover  factor  is  added  based  on  the  cropping-management 
factor  C  used  in  the  USLE  (Wischmeier  and  Smith,  1978).  The  factor  C  is 
defined  as  the  ratio  of  soil  loss  from  land  cropped  under  specified 
conditions  to  the  corresponding  loss  from  tilted,  continuous  fallow. 


Wischmeier  (1972)  presented  a  method  and  the  necessary  graphs  for 

determining  C  for  situations  were  data  are  not  readily  available.  In  this 

method  C  is  treated  as  the  product  of  three  subfactors  depending  on  (i) 

canopy  cover,  (ii)  mulch  and  ground  cover,  and  (iii)  residual  effects  of 

the  land  use,  Cjjj,  respectively.  Approximating  the  first  two  subfactors 

as  linear  functions  of  the  densities  D  and  D  introduced  in  section  2.1, 

c  g  * 

the  factor  C  can  be  expressed  as 

C  =  0.2Cin((l-Dc)(l-Dg)  +  HcDc(l-Dg)] 

where  is  a  function  of  the  average  fall  height  of  drops  from  the  canopy 

cover.  Introducing  C  in  Eq.  59,  neglecting  the  term  depending  on  H  for 

simplicity,  and  incorporating  the  coefficient  0.20^  into  a^,  yields 

D  =  a  I*  (1-D  ) (1-D  ) .  (60) 

r  r  n  c  g 

Ponded  water  deeper  than  a  critical  depth  cushions  the  impact  of  raindrops 
and  also  diminishes  the  erosion  caused  by  the  dissipation  of  impact  energy. 
Mutchler  and  Young  (1975)  suggested  that  a  water  depth  of  more  than  three 
times  the  median  drop  size  essentially  eliminated  detachment  by  raindrop 
impact.  Consequently,  detachment  can  take  place  only  if  the  raindrops  can 
penetrate  through  the  water  depth,  h,  and  the  thickness  of  detached  soil, 
e,  accumulated  from  previous  events.  Laws  and  Parsons  (1943)  found  that 
the  median  drop  size,  expressed  in  millimeters,  is  related  to  rainfall 
intensity  as 


°50  =  2.23  I 


0.182 


(61) 


Using  this  formula,  the  effect  of  water  pondage  is  described  by  modifying 
Eq.  60  as 


Dr  =  ar  *n  U'^)  (l"Dg)  I  l-(h  +  e)/3D5Q],  if 

(h  +  e)  g  3D50,  (62) 

and  Dr  =  0  otherwise.  A  similar  expression  has  been  proposed  by  Li  (1979). 
The  amount  of  soil  detached  by  raindrop  impact  during  a  time  step  At,  and 
to  be  added  to  the  current  storage  of  detached  soil  is,  therefore, 


Aerk  =  arIJ(l-Dc)(l-D  )(l-(h+e)/3D50Jpk  At,  (h+e)  S  3D5Q,  (63a) 

Aerk  =  0,  (h+e)  >  3D5q,  (63b) 

where  pk  is  the  percentage  of  sediment  material  in  the  kth-size  fraction. 

2. 4. 3. 2  Soil  Detachment  by  Flow:  Erosion  by  overland  flows  usually 

occurs  in  many  single  rills.  However,  for  modeling  applications  rill 

erosion  is  assumed  to  be  uniformly  distributed  overland,  and  erosion  by 

concentrated  flows  is  restricted  to  channels.  Erosion  by  flow  potentially 

occurs  if  the  sediment  load  carried  by  the  incoming  flow  is  less  than  the 

transport  capacity  of  the  flow.  If  the  sediment  load  is  greater, 

deposition  occurs.  Erosion  by  flow  is  assumed  to  occur  at  capacity  rate  if 

no  sediment  is  present  in  the  flow.  But  if  the  transport  capacity  of  the 

flow  is  partially  filled  a  corresponding  depletion  of  the  detached  soil 

available  for  transport  is  computed.  If  the  transport  capacity  is  less 

than  the  available  detached  soil,  transport  is  the  limiting  factor  and  no 

erosion  by  flow  takes  place.  If  the  available  detached  soil  is  less  than 

the  transport  capacity,  additional  soil  is  detached  by  the  flow  to 

compensate  for  the  insufficient  supply.  These  concepts  are  implemented  as 

follows.  The  amount  of  sediment  that  the  flow  can  potentially  carry 

downstream  during  the  time  increment  At  is  transformed  into  an  equivalent 
p 

thickness  Ae  ,  as  explained  in  the  following  section.  If  this  thickness  is 

less  than  the  layer  of  available  detached  soil,  e  (Fig.  8a),  no  detachment 

p 

by  flow  occurs.  If  Ae  >  e,  the  available  detached  soil  is  not  enough  for 

transport  and  the  flow  can  potentially  detach  the  additional  amount  e  - 
p 

Ae  .  However,  the  resistance  of  the  soil  to  the  erosive  forces  of  the  flow 
depends  on  the  soil  properties  and  structure,  as  well  as  on  the  condition 
of  the  land  surface  (Olson  and  Wischmeier,  1963;  Kemper,  1966;  Wischmeier 
and  Mannering,  1969;  Grissinger,  1980).  Consequently,  the  potential 
thickness  of  detachment  by  flow  is  modified  by  a  flow  detachment 
coefficient,  a*,  yielding  the  following  amount  of  soil  detached  by  flow: 


(64) 


**fk  =  *af  (AeP  '  e)  Pk* 

where  ranges  from  zero  to  1.0  depending  on  the  soil  erodibility.  The 
detachment  coefficients  in  Eqs .  63  and  64  are  optimization  parameters  that 
are  calibrated  by  fitting  the  sediment  discharge  rates  to  observed  data. 


2.4.4  Numerical  Procedure  for  Sediment  Routing 

The  equations  of  sediment  supply  and  sediment  mass  conservation  are 
coupled  to  the  water  routing  scheme  using  the  procedure  described  below  to 
track  the  evolution  of  the  sediment  load,  and  to  compute  the  aggradation 
and  degradation  of  the  land  surface  and  stream  channels. 

Eq.  48  is  a  linear  hyperbolic  equation  that  may  be  solved  by  the 
method  of  characteristics.  The  steps  given  below  follow  those  presented  in 
section  2-3.2.  The  total  differential  of  Agk(x,t), 
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(65) 


and  Eq.  48  form  a  system  of  two  equations  in  the  two  unknown  partial 
derivatives  of  Agk>  that  is 

"l  V 
dt  dx 

Since  these  derivatives  do  not  exist  along  the  characteristic  paths,  the 
determinant  of  the  coefficient  matrix  vanishes  at  points  on  the 
characteristics.  This  condition  yields  the  characteristic  equation 
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(66) 


Tt  =  V  *  Q'A 


(67) 


in  agreement  with  the  assumption  that  the  sediment  moves  with  the  same 
velocity  of  the  flow.  One  of  the  invariants  of  this  solution  also  gives 
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dt 


=  g£k  +  gdk 


(68) 
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or,  by  integrating  Eq.  68 


sk 


(x’t)  =  Ask(xo’to)  +  {  (g 


£k 


«dk>«- 


(69) 


where  Ask  (xo>to)  is  an  initial  value  of  Applying  this  equation  to 
two  points  E(x^_j,tj_j)  and  F(x^,t^)  on  a  characteristic  path  (Fig.  9)  and 
expressing  the  result  in  discrete  form,  gives 


^Ask\,j  ^Ask^i-l,j-l  + 


“Wi.j  * 


(70) 


in  which  g^  and  g^  are  assumed  uniformly  distribured  in  the  interval  x^ 

p 

=  x  5  x^.  The  potential  transport  capacity  of  the  flow,  C^,  is  obtained 
from  the  formulas  in  section  2.4.2  using  the  average  of  the  flow 
characteristics  computed  at  E  and  F,  and  the  sediment  properties  at  point 

p  - 

E.  Substituting  (C^  A)  for  the  right-hand-side  of  Eq.  70  and  solving  for 

^dk^.j  yields 


(gdk^i,j  "  At.  fCkA  *  ^Ask^i-1 ,  j-1  ^  '  ^Ak^i, 


(71) 


p 

where  the  overbar  denotes  the  average  flow  area,  and  g^  is  the  potential 

change  in  detached  soil  storage  caused  by  either  erosion  or  deposition, 
p 

Deposition:  If  <  0,  the  potential  carrying  capacity  of  the  flow  is 

less  than  the  sediment  present  in  the  flow.  Consequently,  the  sediment 
load  transport  rate  is 


Q  (V,  j  =  Q  <  .  (72) 

P 

and  the  excess  load,  -  g^i  is  deposited  on  the  bed.  However,  whether  this 
amount  will  reach  the  bed  during  the  time  step  At  depends  on  this  being  not 
less  than  the  average  time  for  the  sediment  particles  to  settle  to  the 
channel  bed.  Data  by  Jobson  and  Sayre  (1970)  and  by  Lean  (1971)  indicate 
that  this  settling  time  may  be  computed  using  the  particle  fall  velocities 
in  the  quiescent  fluid.  Therefore,  the  actual  deposition  of  a  size 
fraction  during  the  interval  At  is  calculated  from 
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or 


«Vi,j  *  ■«**'. 


‘Vi.j  *  'KWij  •  lf  *  <  '• 


where  P  =  w^  At/h.  The  new  bed  elevation  is 


1 


n  <V. 


Z.  .  =  Z.  .  .  +  -  l  . 

i.J  i»J"l  p  k=1  1-* 


ill 


(73a) 

(73b) 


(74) 


p 

Erosion:  If  >  0  the  potential  carrying  capacity  exceeds  the  amount  of 

material  in  transport  and,  therefore,  the  flow  will  entrain  additional 
material.  In  this  case,  one  of  two  possible  events  may  occur  depending  on 

p 

whether  g^  is  equal  to  or  greater  than  the  thickness,  (e^^  _j,  of  the 

available  detached  soil.  Let 


(ek}i,j  =  <’ek)i, j-1  +  ^rk^i.j-l  ' 


(75) 


denote  the  depth  of  detached  soil  in  storage  at  the  end  of  the  time  step 

At . ,  in  which 
J 


K>i,j  =  <4\,j  vf 


(76) 


is  the  depth  of  potential  erosion  by  flow. 

(i)  If  (e.  ) .  .  £  0  the  available  detached  soil  is  enough  to  supply 

k  i  ,  j 

the  sediment  entrained  by  the  flow.  In  this  case  no  erosion  (of  undetached 
soil)  occurs,  the  transport  rate  at  the  end  of  the  time  step  At^  is  equal 
to  the  carrying  capacity  (Eq.  72),  and  (Z,  ) .  .  =  (Z.).  .  . 

(ii)  If  (e,  ) .  .  <  0  the  availability  of  detached  soil  is  less  than 

k  i  ,  j 

the  potential  entrainment,  and  additional  soil  is  detached  by  the  flow. 
The  depth  of  erosion  by  flow  is 


(Ae,,).  . 
fk'i.j 


=  -V'k’i.j 


(77) 
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from  which 


'Wi.j =  + 


(78) 


From  Eqns.  70  and  78,  the  concentration  of  the  updated  load  is 


(A  ,).  .  At. 

(CJ;  4  =  +  -T1  (8.J,- 


*  ij 


+  r  l  ( e -  )  -  +  (Ae  )  .  +  (AefJ;  J. 


k  i ,  j - 1  rk  i , j 


(79) 


Finally,  the  new  elevation  of  the  eroded  bed  is 


n  ^efk^i  i 
Z.  .  =  Z.  .  ,  +  I  - 

k=l  l-\ 


(80) 


Substituting  Eq.  72  or  79  into  Eq.  44  gives  the  concentration  of  the  total 
sediment  load.  The  total  concentrations  existing  on  all  the 
characteristics  paths  reaching  a  certain  location  define  the  outflow 
sedimentgraph  at  that  location. 

In  some  instances,  the  time  step  size  selected  for  water  routing  may 
yield  a  cluster  of  small  space  increments  (Fig.  8b).  Repeating  the 
sediment  calculations  for  each  of  these  increments  may  not  be  necessary 
when  the  sediment  routing  parameters  are  only  required  at  points  spaced 
farther  apart.  In  these  cases  it  is  more  efficient  to  combine  a  number  of 
space  and  time  steps  used  for  water  routing,  into  larger  steps  for  sediment 
routing.  To  this  end,  the  program  contains  an  option  to  compute  the 
sediment  routing  parameters  only  at  points  satisfying  the  conditions 


n 

x  =  x  +  I  Ax . 
n  m  i 

i=m 


such  that 


n 

I  Ax.  S  GDX, 
i  =m 


t 

n 


t  + 
m 


n 

l  At. 
J 


where  GDX  is  a  user  supplied  parameter  (Fig.  8b). 


(81a) 


(81b) 
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APPLICATIONS 


To  illustrate  the  relative  accuracy  of  the  present  model,  it  is 
applied  to  several  examples  reported  in  the  literature.  These  are  the 
hypothetical  three-plane  cascade  discussed  by  Kibler  and  Woolhiser  (1970); 
the  experiments  performed  by  Iwagaki  (1955)  involving  unsteady,  open- 
channel  flow  with  lateral  inflow;  the  runoff  studies  carried  out  by  Schaake 
(1970)  in  a  small  urban  catchment;  and  two  agricultural  catchments  reported 
by  Burford  and  Clark  (1973).  The  PSF  and  ASF  schemes  are  compared  on  the 
three  plane  cascade  example.  The  results  from  the  ASF  shock-fitting  scheme 
are  also  compared  with  those  obtained  with  the  implicit  finite-difference 
(FD)  scheme  presented  by  Li  et  al.  (1975a). 

The  computational  parameters  used  in  each  example  are  given  in  Table 
1.  The  kinematic  wave  parameter  n,  Eq.  26,  is  kept  fixed  at  1.50  (Kibler 
and  Woolhiser,  1970;  Singh,  1975).  With  n  fixed,  the  calculations  only 
required  characterization  of  the  parameter  a.  The  values  of  a,  estimated 
by  Kibler  and  Woolhiser  (1970),  are  used  in  the  three-plane  cascade 
computations.  In  all  the  other  examples,  this  parameter  is  computed,  using 
the  relationship  proposed  by  Singh  (1975) 

a  =  Cj  +  C2  (SQ)*  ,  (82) 

where  C^  and  C^  are  constants  to  be  optimized  in  each  case.  This 
relationship  was  chosen  only  for  its  simplicity.  However,  any  other 
expression  could  have  been  used,  since  the  computational  efficiency  of  the 
schemes  is  not  contingent  on  the  manner  or  is  estimated. 

3.1  HYPOTHETICAL  KINEMATIC  CASCADE 

Kibler  and  Woolhiser  (1970)  used  a  three-plane  kinematic  cascade  to 
illustrate  the  formation  of  kinematic  shock  waves.  The  planes  are  400  ft. 
long  with  slopes  0.04,  0.01,  and  0.0025,  respectively.  The  lateral  inflow 
is  a  rainfall  pulse  with  an  intensity  of  0.75  inch/hr  and  a  duration  of  30 
min.  The  necessary  condition  (Eq.  37)  is  satisfied  at  the  unions  of  the 
three  planes.  Thus,  shock- forming  regions  exist  in  the  rising  limbs  of  the 
upstream  inflows  to  the  second  and  third  planes. 

The  characteristics  and  shock  paths  crossing  the  x-t  domain  of  the 
cascade  are  presented  in  Fig.  10.  The  shock  paths  were  computed  using  the 
PSF  scheme.  Because  there  is  no  upstream  inflow  to  the  first  plane,  the 
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outflow  from  this  plane  is  a  smooth  hydrograph  generated  by  the 
characteristic  waves  emanating  from  dry  ground  conditions.  This  hydrograph 
is  used  as  inflow  to  the  second  plane.  The  rising  limb  satisfies  the  shock 
forming  conditions  (Eq.  36),  and  so,  small  shocks  are  introduced  by 
discretizing  this  portion  of  the  hydrograph.  These  shocks  are  projected  to 
the  downstream  boundary  along  with  the  characteristics  emanating  from  the 
initial  dry  ground  condition.  There  they  define  an  outflow  hydrograph 
containing  a  smooth  rising  limb  formed  by  the  characteristic  waves, 
followed  by  a  piecewise  continuous  part  formed  by  the  shock  waves.  To 
reduce  the  number  of  shocks  to  be  traced  in  this  example,  shocks  arriving 
within  the  same  time  step  are  assumed  to  merge  and  proceed  as  a  single 
shock  to  the  next  plane.  The  same  procedure  is  applied  in  routing  over  the 
third  plane.  The  hydrograph  at  the  cascade  outlet  contains  all  the  shocks 
crossing  the  last  two  planes.  The  shocks  slow  down  as  they  enter  the  third 
plane,  reflecting  the  decrease  in  the  value  of  a.  The  complete  outflow 
hydrograph  is  compared  in  Fig.  1 1  to  the  analytical  solution  of  Kibler  and 
Woolhiser  (1970).  Their  solution  displays  the  two  shocks  emanating  from 
the  unions  of  the  planes.  The  two  hydrographs  are  quite  close,  which  shows 
that  the  approximation  introduced  in  formulating  the  PSF  solution  does  not 
detract  from  its  accuracy. 

The  same  hydrograph  calculation  was  carried  out  using  the  ASF  and  FD 
schemes.  The  ASF  solution  is  plotted  in  Fig.  11.  It  exhibits  a  smooth 
rising  limb  showing  the  overall  effect  of  the  shock  waves,  and  is  in  good 
agreement  with  the  two  other  solutions.  Moreover,  it  is  twice  as  fast  as 
the  PSF  scheme  (Table  1).  It  is  thus  evident  that  the  ASF  scheme,  in 
addition  to  being  easier  to  work  with,  gives  results  as  accurate  as  the 
analytical  solution.  For  these  reasons  the  PSF  technique  was  abandoned  in 
favor  of  the  more  efficient  ASF  solution. 

The  hydrograph  computed,  using  the  FD  scheme  with  the  same  values  of 
a,  is  compared  in  Fig.  12  with  the  ASF  solution.  Although  the  hydrographs 
agree  in  their  predicted  yields,  the  FD  solution  exhibits  a  pronounced 
delay  and  reduction  of  the  peak.  In  an  attempt  to  improve  this  solution, 
the  FD  calculation  was  repeated  using  a  new  optimized  value  of  the 
kinematic  parameter.  The  result  shows  improvement  in  the  peak  estimate  but 
an  overall  deterioration  of  the  hydrograph  shape  (Fig.  12).  This  inability 
of  the  FD  scheme  to  reproduce  the  analytical  solution  is  further  discussed 
in  the  context  of  the  next  example. 
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Fig.  12  Comparison  of  finite-difference  results  with  approximate 
shock-fitting  solution 
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UNSTEADY  CHANNEL  FLOW 


Iwagaki  (1955)  reported  a  set  of  laboratory  experiments  of  unsteady 
flow  with  lateral  inflow  in  a  smooth  open  channel  7.3  ft  long.  The  channel 
was  divided  into  three  sections  of  equal  lengths  and  different  slopes. 
From  the  upstream  end,  the  slopes  were  0.020,  0.015,  and  0.010.  The 
lateral  inflow  was  adjusted  so  that  each  section  would  receive  a  constant 
rate  of  inflow,  with  the  upper,  middle,  and  lower  receiving  0.0425,  0.0251, 
and  0.0315  inch/sec.  The  durations  of  lateral  inflow  used  in  the 
experiments  were  10,  20,  and  30  sec.  Under  these  flow  conditions,  Eq.  37 
predicts  a  shock  forming  zone  at  the  union  of  the  upper  and  middle 
sections . 

Figs.  13,  14,  and  15  show  that  the  hydrographs  computed  with  the  ASF 
scheme  are  in  good  agreements  with  the  partial-equilibrium  hydrographs, 
measured  by  Iwagaki.  The  kinematic  parameter  a  was  adjusted  by  fitting  the 
data  plotted  in  Fig.  13.  Because  the  slopes,  geometries,  and  hydraulic 
roughness  were  the  same  in  all  three  experiments,  the  same  value  of  a  was 
used  in  the  calculation  of  the  other  two  hydrographs.  These  solutions 
correctly  simulated  the  shock  formation  and  its  arrival  time  at  the  outlet. 
This  is  clearly  depicted  in  Fig.  15,  where  the  hydrograph  peaks  immediately 
after  the  lateral  inflow  terminates  and  falls  continuously,  until  a  sudden 
rise  occurs  when  the  shock  arrives.  Fig.  14  shows  the  shock  arriving 
shortly  after  the  first  peak,  whereas  in  Fig.  13  the  shock  arrives  well 
ahead  of  the  peak. 

Two  different  hydrographs  computed  with  the  FD  scheme  are  given  in  the 
same  figures.  The  hydrographs  computed  with  the  values  of  a  used  ii.  the 
ASF  solution  do  not  agree  very  well  with  the  measurements  and,  in  addition, 
they  do  not  reproduce  the  aforementioned  shock  effect.  A  second  solution 
was  obtained  by  recalibrating  the  kinematic  parameter  but  the  new 
hydrographs  do  not  exhibit  any  better  agreement.  In  the  present  example, 
the  shock  formation  is  an  intrinsic  part  of  the  solution  to  the  kinematic 
wave  approximation.  However,  because  of  its  smoothing  effect,  the  FD 
scheme  is  unable  to  reproduce  this  important  aspect  of  the  analytical 
solution.  In  addition,  this  scheme  required  between  50  and  100  percent 
more  computing  time  that  the  ASF  solution  (Table  1). 
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Fig. 


13  Simulation  of  Iwagaki's  partial-equilibrium  hydrograph  for  lateral 
inflow  duration  of  30  sec 


Fig.  14  Simulation  of  Iwagaki's  pa rt ial-equi 1 ibri um  hydrograph  for  lateral 
inflow  duration  of  20  sec 


ig.  15  Simulation  of  Iwagaki's  partial-equilibrium  hydrograph  for  lateral  inflow  duration  of  10  sec 
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URBAN  CATCHMENT 


This  example  is  used  to  illustrate  that  the  shock-fitting  technique 
presented  in  this  report  is  applicable  to  a  combination  of  interconnected 
overland  and  channel  segments.  The  case  considered  herein  is  the  small 
urban  catchment,  reported  by  Schaake  (1970).  In  this  publication  Schaake 
describes  an  event  labeled  3  SPL1  on  a  0.39  acre  impervious  parking  lot 
labeled  SPL1.  He  represented  the  catchment  by  a  number  of  interconnected 
overland  segments  and  V  shaped  channels  (Fig.  16).  The  overland  segments 
vary  from  20  to  36  ft  in  length,  and  their  slopes  range  from  0.0167  to 
0.019.  The  V-shaped  channels  have  lengths  varying  from  50  to  165  ft,  and 
slopes  ranging  from  0.0148  to  0.0213.  The  channel-side  slopes  were  all 
1:11.3.  The  same  geometric  representation  is  used  in  the  simulations 
reported  herein. 

The  runoff  hydrographs  measured  at  the  outlet  of  the  catchment  and  the 
hydrographs  computed  with  the  ASF  and  FD  schemes  are  shown  in  Fig.  17.  The 
value  of  a  was  obtained  by  fitting  the  measured  data.  The  result  obtained 
with  the  ASF  scheme  overpredicts  the  first  peak  but  agrees  very  well  with 
the  rest  of  the  data.  The  FD  calculation  obtained  with  the  same  a  does  not 
predict  the  second  peak  well.  In  an  attempt  to  improve  the  second  peak,  a 
second  calculation  using  the  FD  scheme  was  made  by  recalibrating  a.  This 
improved  the  second  peak  but  did  not  help  the  rest  of  the  hydrograph  (Fig. 

17) .  In  general,  the  hydrographs  predicted  by  the  FD  scheme  do  not  exhibit 
the  significant  deviations  shown  in  the  previous  calculations.  This  is 
because,  in  this  example  the  lateral  inflow  rate  changes  rather  smoothly, 
thus  reducing  the  magnitudes  and  effects  of  shocks.  Also,  the  slopes  in 
the  previous  two  examples  exaggerate  the  development  of  shocks,  whereas  in 
this  example  the  slopes  are  nearly  the  same,  thus  leading  to  much  smaller 
shocks.  Nevertheless,  these  solutions  take  twice  as  much  computing  time  as 
the  ASF  calculation  (Table  1). 

3.4  AGRICULTURAL  CATCHMENTS 

The  data  used  in  these  tests  were  obtained  from  the  USDA  experimental 
catchments  W-5  southwest  of  Holly  Springs,  MS,  and  R-5  near  Chickasha,  OK 
(Burford  and  Clark,  1973).  Catchment  W-5  drains  a  1.76  mile2  area  (Fig. 

18) ,  with  a  good  mixturp  of  cultivated  land,  timber,  pasture,  and  idle 
land.  Catchment  R-5  has  an  area  of  23.7  acres,  and  is  range  land  with  an 
excellent  native  grass  cover  (Fig.  19).  These  catchments  were  chosen 
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because  of  diversity  and  availability  of  most  of  the  required  data. 
Sediment  records  were  available  for  catchment  R-5,  but  the  sediment  yield 
was  so  small  that  no  comparisons  were  made  between  measured  and  computed 
sediment  discharges. 

The  model  was  calibrated  using  the  event  of  February  21,  1971,  on 
catchment  W-5,  and  the  event  of  May  6,  1969,  on  catchment  R-5. 

Infiltration  parameters  for  catchment  W-5  were  estimated  from  information 
reported  by  Smith  and  Parlange  (1978)  for  Colby  swelling  type  soils, 
because  very  little  infiltration  data  were  available.  The  infiltration 
parameters  for  catchment  R-5  were  obtained  from  an  average  infiltration 
curve  obtained  from  a  large  number  of  inf iltrometer  runs  conducted  in  the 
fall  of  1977.  The  calibrations  were  carried  out  by  adjusting  first  the 
flow  resistance  parameters  (Eq.  82)  to  match  the  hydrographs,  and  then  the 
sediment-model  parameters  (Eqns .  63  and  64)  were  adjusted  to  fit  the 

sedimentgraphs .  These  same  parameters  were  used  in  simulating  all  the 
remaining  events.  The  results  of  these  tests  were  reported  in  an  earlier 
paper  by  Alonso  et  al.  (1978). 

Examples  of  the  comparison  between  the  simulated  and  the  measured 
hydrographs  and  sedimentgraphs  are  shown  in  Figs.  20,  21,  and  22.  A  total 
of  nine  events  on  catchment  W-5  and  two  on  catchment  R-5  were  simulated. 
The  agreement  between  the  shapes  of  measured  and  simulated  events  is 
satisfactory.  Comparisons  between  measured  and  computed  water  yield,  peak 
runoff  rates,  sediment  yield,  and  peak  sediment  discharge  are  given  in 
Figs.  23  and  24.  These  plots  show  that  the  model  estimates  both  yields  and 
peaks  within  a  range  of  about  ±40  percent  of  measured  values.  The  limited 
number  of  events  used  in  this  study  precluded  an  estimation  of  the 
confidence  level  of  the  range  of  variability.  The  results  presented  in 
Figs.  20  through  24  indicate  that  simulations  of  different  size  events 
using  only  one  set  of  parameters  for  each  catchment,  were  satisfactory. 
This  suggests  that  the  model  could  be  used  to  predict  the  response  of  a 
catchment  to  different  management  practices,  if  the  model  parameters 
associated  with  each  practice  could  be  accurately  estimated.  Also,  the 
above  results  indicate  that  the  model  could  be  readily  transferred  to 
ungaged  catchments,  if  the  model  parameters  were  properly  regionalized. 
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Fig.  21  Sedimentgraph  from  catchment  W-S  for  the  February  21,  197]  event 
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24  Comparison  of  measured  and  computed  sediment  yields  and  peaks 
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CONCLUSIONS  AND  RECOMMENDATIONS 


4 

4.1  CONCLUSIONS 

1.  A  numerical  model  for  routing  water  and  sediment  on  small 
catchments  has  been  developed.  The  model  accepts  any  single  rainfall 
hyetograph  and  produces  runoff  and  sediment  hydrographs  for  the  modeled 
catchment. 

2.  The  model  is  developed  on  a  general  basis  so  that  it  may  be 
applied  to  any  agricultural  catchment  by  changing  only  the  input  data.  The 
approximate  range  of  basins  over  which  the  model  is  applicable  is  from  a 
few  acres  to  about  5  square  miles. 

3.  The  model  is  based  on  the  physical  processes  governing  the 
mechanics  of  water  and  sediment  movement  and  requires  the  calibration  of 
four  parameters. 

4.  The  model  can  be  used  to  simulate  the  effect  of  different  land 
uses  on  the  water  and  sediment  yields  from  the  modeled  catchment. 

5.  The  model  predicts  the  surface  component  only.  It  does  not 
presently  predict  subsurface  and  groundwater  movement. 

6.  The  applicability  of  the  model  is  restricted  to  streams  where  the 
channel  geometry  does  not  change  significantly  during  a  storm  event,  and  in 
which  the  kinematic-wave  approximation  for  flow  routing  is  valid. 

7.  The  present  model  simulates  single  storm  events;  the  user  must 
estimate  the  initial  conditions  for  the  storm.  The  model  can  still  be 
applied  to  a  sequence  of  rainfall  events  if  the  user  can  make  satisfactory 
estimates  of  the  initial  soil  moisture  conditions.  In  this  case  the  model 
can  be  used  to  predict  a  sequence  of  surface  runoff  and  sediment  transport 
events . 

8.  The  model  has  been  validated  with  several  data  sets  including  data 
from  the  natural  catchments  W-5  in  northern  Mississippi,  and  R-5  near 
Chickasha,  Oklahoma.  The  shape  of  water  and  sediment  hydrograph  and  total 
water  and  sediment  yields  of  a  number  of  events  were  satisfactorily 
simulated. 

4.2  RECOMMENDATIONS 

1.  It  is  recommended  that  the  model  be  put  to  work  on  real  systems. 
The  evaluation  and  continuous  updating  of  the  model  are  essential  to  its 
credibility  and  effectiveness. 
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2.  It  is  recommended  that  the  model  be  further  developed,  or 
modified,  to  permit  continuous  simulation  over  long  time  spans  (20  to  50 
years).  This  is  essential  in  evaluating  the  long  term  response  of  a 
catchment,  or  stream  network,  which  is  dependent  not  only  on  the  history  of 
management  practices,  but  also  on  the  sequence  of  storm  events. 

3.  It  is  recommended  that  the  model  be  further  developed  to  track  the 
channel  geometry  of  streams  that  become  unstable  due  to  bank  erosion  and 
deposition. 

4.  It  is  recommended  that  data  gathering  efforts  be  continued  to 
provide  an  adequate  base  for  further  model  development  and  validation. 
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ADDENDUM  1.  DESCRIPTION  OF  DATA  INPUT  AND  COMPUTER  PROGRAM 


The  structure  of  the  program  SEDLAB,  a  software  system  developed  for 
the  simulation  of  water  and  sediment  movements  in  agricultural  catchments 
is  described.  The  numerical  schemes  on  which  the  program  is  based  have 
been  described  in  detail  in  Part  2.  A  description  of  the  data  input  and 
important  variables  used  in  the  program  is  given  in  the  following  sections 
of  this  addendum.  The  last  section  presents  a  complete  list  of  the 
program. 

The  system  SEDLAB  consists  of  a  main  program  and  several  subroutines. 
The  sequence  of  operations  performed  by  the  main  program  is  shown  in  Fig. 
1.1.  The  main  program  inputs  the  required  data  to  the  system,  calls  the 
subroutines  according  to  the  computation  scheme,  and  prints  out  the 
calculated  results.  Subroutine  INTCPN  computes  the  evaporation  and 
interception  losses,  and  determines  the  net  rainfall  rate.  Subroutine 
INFLTN  computes  the  ponding  time,  and  the  infiltration  losses  in  a  segment 
for  one  time  step.  Subroutine  WROUT  routes  flow  through  a  segment  for  one 
time  step  and  it  calls  in  turn  subroutine  SROUT.  This  subroutine  performs 
sediment  routing  through  a  segment  for  one  time  step.  Subroutines  LAURSN, 
SETVEL,  SHIELD,  YALIN,  and  YANG  are  called  by  SROUT  to  compute  potential 
carrying  capacities  and  sediment  transport  parameters.  The  program  has 
been  designed  to  process  the  entire  network  of  segments  for  one  time  step 
according  to  the  computational  sequence  described  in  section  1.2  of  this 
addendum.  Once  the  entire  network  has  been  processed  the  simulation  neves 
to  the  next  time  step.  This  sequence  is  repeated  until  the  entire 
simulation  period  is  computed. 

The  computer  code  requires  49,710  words  on  a  Mod  Comp  Classic  computer 
system.  This  is  a  16-bit  machine  that  uses  two  words  for  each  single 
precision  variable.  The  execution  time  for  the  event  of  February  21,  1971, 
on  catchm>  v-5  is  90  seconds. 

A!  1  DATA  INPUT 

'uta  inquired  to  run  program  SEDLAB  includes: 

«•  !(•  i  i  c.t,  length,  bed  slope,  bed  elevation,  and  wetted 

•  .  i.  ,  •  i-  *  i  »■  1  at  i  ons  . 

.  •  .  if.  (  .  and  ground  cover  density,  interception 
••  !  ground  (overs,  and  ratio  of  evaporating 
,  -  I  f  :  g  1 1  und  r uve r . 


Soil  data:  soil  saturated  hydraulic  conductivity,  soil  sorptivity, 
specific  gravity  and  si2e  distribution  of  bed  material. 

Flow  and  sediment  routing  parameters:  constants  describing  flow 

resistance,  parameters  describing  soil  detachment  by  raindrop  impact  and 
surface  runoff,  maximum  penetration  depth  of  raindrop  impact,  and 
computational  sequence. 

Storm  characteristics:  rainfall  intensity,  mean  evaporation  rate,  and 
initial  interception  storage  content. 

A  detailed  description  of  the  data  input  is  given  below: 


Card  FORTRAN 

No.  Variable 

Descript  ion 

Units 

Format 

1  TITLE 

Alphanumerical  identification  of 

simulation  run. 

- 

20A4 

2  AREA 

Drainage  area  of  catchment 

acres 

F10.4 

NOV 

Total  number  of  overland  segments 

- 

14 

NCH 

Total  number  of  channel  segments 

with  negligible  infiltration 

"" 

14 

NCI 

Total  number  of  channel  segments 

with  significant  infiltration 

— 

14 

NSTHM 

Total  number  of  storm  events 

simulated  in  the  run 

14 

3  (This  card  must  he  repeated  lor  each  overland 

segment ) 

SEG 

Number  ident living  overland  segment 

Segments  are  numbered  sequentially 

from  1  to  NOV. 

14 

OVA 

Area  of  overland  segment 

ft* 

F12.3 

SEEN 

Slope  length  of  overland  segment. 

This  length  is  computed  by  dividing 

OVA  by  the  length  of  the  receiving 

channel  reach. 

ft 

F10.4 

SLOPE 

Slope  of  overland  segment. 

ft/ft 

F10.4 

CPER 

Coefficient  in  the  wetted  perimeter 

versus  flow  area  relation  (the 

default  value  is  one) 

F10.4 

I  .83 


EPER 


F10.4 


Exponent  in  the  wetted  perimeter 
versus  flow  area  relation  (the 

default  value  is  zero) _ 

(This  card  must  be  repeated  for  each  channel  segment) 

SEG  Number  identifying  channel  segments.  -  14 

These  segments  are  numbered  sequentially 
from  (NOV+1)  to  (NOV+NCH)  if  there  are 
no  channel  reaches  with  significant 
infiltration.  Otherwise,  the  channel 
segments  are  numbered  from  (NOV+1)  to 
(NOV+NCH+NCI). 


SLEN 

SLOPE 

CPER 

EPER 

Length  of  channel  segment 

Bed  slope  of  channel  segment 

Coefficient  in  the  wetter  perimeter 

versus  flow  area  relation. 

Exponent  in  the  wetted  perimeter 

versus  flow  area  relation. 

ft 

ft/ft 

F10.4 

F10.4 

F10.4 

F10.4 

VOG 

Interception  storage  capacity  for  a 

typical  ground  cover. 

inches 

F10.4 

SRG 

Ratio  of  evaporating  surface  to  the 

horizontal  projected  area  of  typical 

ground  cover 

it1 lit2 

F10.4 

VOR 

Ratio  of  the  interception  storage 

capacity  of  a  typical  canopy  cover  to 

that  of  a  typical  ground  cover. 

F10.4 

HLR 

Average  height  of  ground  cover  in 

channels . 

ft 

F10.4 

NSOIL 

Number  of  representati ve  sediment 

size  fractions  used  in  the  simulation. 

- 

14 

SPGR 

Specific  gravity  of  sediment 

- 

F10.4 

AIM 

Coefficient  of  soil  detachment  by 

- 

F10.8 

raindrop  impact  (a^).  User  supplied 
optimization  parameter. 


ADF  Coefficient  of  soil  erosion  by  surface 

flow  (a^).  User  supplied  optimization 
parameter. 


F10.8 


GMAX 

Maximum  penetration  depth  of 

raindrops  (Eq.  61). 

ft 

F10.4 

GDX 

Space  increment  adopted  for  sediment 

routing.  User  supplied  parameter. 

ft 

F10.4 

7 

D50 

Median  size  of  sediment  fractions. 

mm 

10F7.3 

The  program  can  accomodate  up  to 

five  fractions. 

8 

PC 

Percentage  of  sediment  fractions. 

* 

10F7.3 

9 

(This 

card  must  be  repeated  for  each  segment  in 

the  sequences  used 

for  ' 

cards  3  and  4) 

CANO 

Canopy  cover  density  for  the  segment 

ft2/ft2 

F10.4 

GCOV 

Ground  cover  density  for  the  segment 

ft2/ft2 

F10.4 

HYCND 

Saturated  hydraulic  conductivity  for 

the  segment 

in/hr . 

F10.4 

10 

(This 

card  must  be  repeated  for  eai h  segment) 

ISEG 

Index  identifying  the  position  of  the 

segment  in  the  computational  sequence 

(see  following  section) 

14 

I  ARY 

Array  containing  the  storage  locations 

- 

514 

of  the  inflows  and  outflows  computed  lor 

the  segment  (see  following  section) 

11 

TEMP 

Water  temperature 

Celsius 

F10.4 

GAMA 

Specifu  weight  uf  w.itci 

lb/ 1 1  3 

no. 4 

SNU 

Kinemalu  visiosity  <>t  water 

( t  2 /sec 

F10.8 

EXP 

Exponent  in  fq.  2<> 

- 

no. 4 

Cl 

First  eoel  lie  lent  desinhing 

- 

F10.4 

ki  Hemal  n  -wave  f  rut  ion  parameter 

(F.q.  82  1.  User  supplied  optimization 

pa  ramete r 

C2 

Second  coefficient  in  Eq.  82. 

User  supplied  optimization  parameter 

- 

F10.4 

12 

NEED(I)  Vector  representing  the  following 

- 

614 

output  options: 
NEED  (1)  input  data 
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NEED  (2)  Bed  elevation  changes  computed  for 
each  segment  at  the  end  of  the 
simulation  period 

NEED  (3)  water  budget  and  sediment  yield 
NEED  (4)  Computed  infiltration  rates.  Measured 

and  computed  hydrographs  and  sedimentgraphs 
NEED  (5)  Plots  of  hyetograph  and  hydrograph 
NEED  (6)  Plots  of  hyetograph  and  sedimentgraph 


When  NEED(I)  >  0,  the  program  prints  selected  output; 
when  NEED(I)  =  0  the  program  does  not  print  selected 
output . _ 


13  (Cards 

13  through  17  must  be  repeated  for  each 

storm  event) 

STORM 

Alphanumerical  identification  of  the 

storm  event 

5A4 

DTM 

Size  of  time  step 

min 

F10.4 

ITMAX 

Number  of  time  steps  in  rainfall 

duration 

14 

ITCOM 

Number  of  time  steps  in  simulation 

period 

14 

EVP 

Mean  evaporation  rate 

in/hr 

F10.4 

VIN 

Initial  interception  storage,  defined 

in/in 

F10.4 

as  the  ratio  of  the  initial  storage 
to  the  intercept  ion  storage  capacity 


14 

(This  card  must  be  repeated  every  eight  segments  until  all 

segments  are  included) 

SORPTY  Sorptivity  parameter  for  each  of  the  in2/hr 

eight  consecutive  segments  (Sz/2, 

Fig.  3) 

the 

8F7.7 

15 

(This  card  must  be  repeated  every  twelve  time 

rainfall  intensities  are  read  in) 

steps  until 

all  the 

DR  Rainfall  intensity  for  each  of  the 

twelve  consecutive  time  steps 

in/hr 

12F6.3 

16 

(This  card  must  be  repeated  every  twelve  time 

steps  are  included) 

steps  until 

all  ITCOM 

QMES  Water  discharge  measured  at  the 

catchment  outlet  for  each  of  the 

twelve  consecutive  time  steps 

ft3/sec 

12F6.2 
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17  (This  card  must  be  repeated  every  twelve  time  steps  until  all  ITCOM 
steps  are  included) 

GMES  Total  sediment  discharge  measured  at  Ibs/sec  12F6.2 

the  catchment  outlet  for  each  of  the 

_ twelve  consecutive  time  steps _ 

A1.2  PREPARATION  OF  INPUT  ARRAYS  ISEG  AND  I ARY 

The  first  step  in  preparing  these  data  is  to  divide  the  catchment  into 
interconnected  overland  and  channel  segments,  each  homogeneous  within 
itself,  and  to  assign  an  identification  number  to  each  segment.  The 
principal  direction  of  flow  is  determined  for  each  overland  segment  from 
the  contour  line  map  of  the  catchment.  The  flow  path  through  the  cascade 
of  segments  is  established  by  following  the  logics  of  gravity  and  flow 
continuity.  The  order  in  which  the  segment  numbers  appear  in  the  flow  path 
defines  the  computational  sequence.  This  procedure  is  illustrated  in  Fig. 
1.2  which  presents  the  segmentation  used  in  simulating  the  catchment  W-S 
described  in  section  3.4.  The  arrows  shown  in  the  figure  denote  the  flow 
direction  in  each  overland  segment.  It  should  be  noticed  that  the 
segmentation  is  restricted  to  no  more  than  two  overland  flow  segments  as 
input  to  a  channel  segment,  and  at  most  two  inflow  channel  segments  to  a 
downstream  channel  segment. 

After  the  segmentation  of  the  catchment  has  been  completed,  two  arrays 
must  be  6et  up  by  the  user.  The  first,  [SEG,  tells  the  program  the 
sequence  for  computing  flow  and  sediment  discharge  for  each  segment.  The 

other,  IARY,  tells  the  program  for  any  segment  where  to  find  previously 

computed  inflows  and  where  to  store  the  computed  outflows.  These  inflows 
and  outflows  are  stored  in  different  columns  of  the  matrices  Q  and  GS 
described  in  the  next  section. 

A  table,  AUX  (Fig.  1.3),  is  used  as  an  aid  to  set  up  the  two  arrays. 
AUX  and  the  matrices  Q  and  GS  have  the  same  number  of  columns.  The  table 
is  constructed  as  follows.  Starting  with  one  of  the  most  upstream  channel 
segments,  its  inflow  segments  are  selected.  Then  their  numbers  are  placed 
in  the  first  available  columns  of  AUX  and  in  separate  rows.  This  will 

usually  mean  an  overland  flow  segment  number  in  row  1,  column  1,  and 

another  in  row  2,  column  2.  The  channel  segment  number  is  then  placed  in 
the  next  available  row  and  the  first  available  column  if  no  further  inflows 
to  the  segment  need  to  be  computed.  An  available  column  is  one  that  does 


1.87 


LEGEND 

-  CATCHMENT  BOUNDARY 

OVERLAND  SEGMENT  ~ 

CHANNEL  SEGMENT 

/  dfS^  \©lx 

-  ^/tThO®/ 

/  y  |i%  / 

TV  at  7  / — / 

scxxy/y  X®  / 


®—  S(ikiX|y 


/ 


H/ 


&'  ©' 

/  VV 

/ 


X  // 


/  @ 


OUTLET 


1.2  Geometric  segmentation  and  flow  path  for  catchment  W-5 


not  have  the  number  of  a  segment  waiting  to  be  combined  in  some  further 
downstream  segment.  This  procedure  is  continued  through  the  flow  path, 
inserting  segment  numbers  in  the  table  until  the  outlet  of  the  catchment  is 
reached.  If  a  junction  with  another  channel  segment  is  reached,  then  what 
has  been  computed  must  be  retained  while  the  area  upstream  of  the  junction 
is  computed. 

To  illustrate  consider  Fig.  1.2.  Channel  segment  29  is  chosen  as  the 
first  most  upstream  channel  segment.  Other  segments  which  could  have  been 
chosen  as  well  to  start  are  30,  33,  34,  35,  or  20.  The  overland  flow 
segments  for  segment  29  are  1  and  2.  These  are  placed  respectively  in  row 
1,  column  1,  and  row  2,  column  2.  When  the  two  overland  flow  segments  are 
combined  to  give  the  channel  outflow  in  29,  they  are  no  longer  needed. 
This  allows  the  outflow  from  29  to  be  t  tered  in  column  1  of  the  next  row. 
Next,  the  outflow  from  channel  segment  30  must  be  computed.  This  requires 
computing  the  overland  flow  segments  3  and  4.  They  are  entered  in  columns 
2  and  3  because  29  must  be  saved  and  is  in  the  column  1.  When  30  is 
computed  it  must  be  saved  also  until  the  inputs  5  and  6  to  segment  31  have 
been  computed.  When  29,  30,  5,  and  6  have  all  been  computed,  then  31  can 
be  computed  and  placed  in  the  first  column  of  the  next  row.  None  of  the 
segments  previously  computed  need  to  be  retained  at  this  point.  When 
segment  32  has  been  computed  it  must  be  retained  until  segment  39  has  been 
computed  which  requires  moving  to  upstream  segments  above  39.  The  channel 
segments  12,  33,  and  38,  and  their  corresponding  lateral  inflows,  are  the 
upstream  starting  points  above  39. 

In  the  given  example,  the  segment  numbers  entered  in  AUX  occupy  five 
columns.  However,  a  finer  segmentation  or  a  denser  drainage  network  would 
have  resulted  in  a  larger  number  of  columns.  The  number  of  rows  of  AUX  is 
always  equal  to  the  total  number  of  segments  because  each  segment  is 
computed  only  once  in  each  time  step. 

When  AUX  is  completed,  ISF.G  and  IARY  are  constructed  from  it  and  then 
AUX  is  discarded.  To  construct  the  column  matrix  ISEG,  the  number  of  the 
segment  in  each  row  of  AUX  is  placed  in  the  corresponding  row  of  ISEG  (Fig. 
1.3).  The  order  of  segment  numbers  in  ISEG  gives  the  computational 
sequence  used  by  the  program  to  route  water  and  sediment  through  the  flow 
path. 
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Fig.  1.3  Assemblage  of  input  arrays  ISEG  and  I ARY  for  catchment  W-5 
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To  construct  IARY,  each  segment  in  the  catchment  is  considered.  The 
information  for  segment  number  n  is  placed  in  the  nth-row  of  IARY.  Columns 
1  and  2  of  IARY  contain  information  about  any  overland  flow  segments  which 
are  input  to  that  segment.  These  columns  will  contain  a  zero  if  there  are 
no  overland  flow  segments,  and  will  contain  the  column  locations  for  the 
overland  flow  segments  in  AUX  if  there  are  (Fig.  1.3).  For  example, 
segment  16  is  an  overland  segment  receiving  no  overland  contributions 
itself.  Thus  columns  1  and  2  of  row  16  in  IARY  have  zeros.  The  lateral 
inflows  to  segment  39  come  from  the  overland  segments  21  and  22  which  are 
found  in  columns  4  and  5  of  AUX.  Thus  columns  1  and  2  of  row  39  in  IARY 
have  the  numbers  4  and  5,  respectively. 

Columns  3  and  4  of  IARY  store  information  about  the  channel  inflows 
for  downstream  channel  segments.  If  the  nth-segment  is  a  channel  segment 
receiving  channel  inflows,  columns  4  and  5  of  the  nth-row  in  IARY  will  have 
the  column  location  of  these  inflow  segments  in  AUX.  For  example,  segment 
12  is  an  overland  flow  segment  and  has  no  channel  inflows,  thus  columns  3 
and  4  of  IARY's  row  12  are  zeros.  Similarly,  segment  29  is  a  channel 
segment  with  no  upstream  channel  inflows,  thus  columns  3  and  4  of  row  29  in 
IARY  are  also  zeros.  Alternatively,  segment  41  has  one  upstream  channel 
inflow  (segment  40)  which  is  found  in  column  1  of  AUX  when  computed. 
Therefore,  column  3  of  row  41  in  IARY  contains  a  1  while  column  4  is  a 
zero.  Segment  39  has  two  channel  inflows  (segments  37  and  38).  When  these 
upstream  segments  are  computed  they  are  entered  respectively  in  columns  2 
and  3  of  AUX.  Therefore,  columns  3  and  4  of  row  39  in  IARY  contain  2  and 
3,  respectively. 

Column  5  of  IARY  contains  the  column  location  in  AUX  of  the  computed 
segment  outflows.  For  example,  the  outflow  from  segment  16  when  computed 
is  entered  in  column  5  of  AUX.  Thus  row  16,  column  5,  of  IARY  contains  a 
5.  Similarly,  the  outflow  of  segment  28  is  found  in  column  3  of  AUX.  Thus 
row  28,  column  5,  of  IARY  contains  a  3. 

The  matrices  ISEG  and  IARY  are  used  sequentially  to  route  water  and 
sediment  through  the  flow  path  for  each  time  step.  For  each  segment  number 
listed  in  ISEG,  the  numbers  in  the  corresponding  row  of  IARY  tell  the 
program  where  the  associated  inflows  are  stored  in  Q  and  GS .  When  all  the 
segments  have  been  worked  through,  the  entire  procedure  is  repeated  for  the 
next  time  step. 
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A1.3  FORTRAN  VARIABLES  NOT  INCLUDED  IN  THE  INPUT  DATA  LIST 


Name 

ACCLN 

ADP 

ALP 

AP 

BEG 

CND 

CONC 

CP 

DFT 

DX 

END 


ERO 

GCD 

GDLAT 

GDUP 

GS 


GTOT 

HYR 


Description 

Gravitational  acceleration 

Thickness  of  detached  soil  at  the  start  of  the 
current  time  step 

Kinematic-wave  friction  coefficient  in  Eq.  26 
Cross-sectional  area  of  flow  at  the  start  of  the 
current  time  step 

Distance  of  space  increment  used  in  sediment 
routing  to  the  upstream  boundary  of  the  segment 
Canopy  cover  density 

Volume  concentration  of  individual  material 
fraction  at  carrying  capacity 
Volume  concentration  of  individual  material 
fraction  at  the  start  of  the  current  time  step 
Medium  size  of  se  .  uit  fraction 
Space  increment  used  in  sediment  routing 
Distance  travelled  by  sediment  characteristic  at 
the  end  of  the  current  time  step,  and  measured 
from  upstream  boundary  of  segment 

Volume  rate  of  soil  detachment  by  raindrop  impact 
Ground  cover  density 

Volume  rate  of  lateral  inflow  of  material  fraction 
to  segment  during  current  time  step 
Volume  rate  of  upstream  inflow  of  material 
fraction  to  segment 

Volume  discharge  of  individual  material  fraction. 
This  is  a  five-column  matrix  used  to  store  the 
sediment  inflow  and  outflow  for  a  segment.  This 
matrix  has  as  many  rows  as  time  steps  in  the 
simulation  period. 

Total  sediment  discharge  at  the  catchment  outlet 
Flow  depth  in  overland  segments,  and  hydraulic 
radius  in  channels 


Units 

ft/sec2 

ft 

ft^/sec 

ft2 

ft 

ft2/ft2 

ft3/ft3 

ft3/ft3 

ft 

ft 


ft 

ft3/ft?  sec 
ft2/ft2 

ft3/ft.  sec 

ft3/sec 


ft3/sec 
lbs/ sec 

ft 


UPON  Number  of  time  steps  until  time  of  ponding  from 
the  start  of  the  simulation  period 

ITYPE  Index  defining  type  of  flow  in  segment.  Set 

ITYPE  equal  to  1  for  overland  flow,  2  for  channel 
flows  with  negligible  infiltration,  and  3  for 
channel  flows  with  significant  infiltration 

KOUT  Column  in  matrix  IARY  containing  storage  location 

of  segment  outflow  computed  at  the  end  of  the 
current  time  step. 

KS  Number  of  time  steps  until  the  time  of  formation 

of  a  characteristic,  or  shock,  wave. 

K!  Number  of  the  current  time  step 

K2  Index  identifying  a  characteristic  wave  emanating 

from  dry  ground  (K2=0),  or  from  upstream  boundary 
(K2>0) 

NO  Number  of  current  space  increment  when  routing 

a  sediment  characteristic 

NSEG  Total  number  of  segments  in  the  catchment 

PERM  Saturated  hydraulic  conductivity  in/hr 

POR  Porosity  of  bed  material  ft3/ft3 

Q  Water  discharge.  This  is  a  five-column  matrix 

used  to  store  the  water  inflow  and  outflow 
computed  for  a  segment.  This  matrix  has  as  many 


rows  as  time  steps  in  the  simulation  period.  ft3/sec 

QE  Water  discharge  at  the  end  of  the  current  sediment 

routing  step.  ft3/sec 

QINF  Rate  of  infiltration  during  the  current  time  step  in/hr 

QLAT  Lateral  inflow  of  water  to  segment  during  the 

current  time  step.  ft3/ft.  sec 

QOUT  Water  discharge  at  the  catchment  outlet  ft3/sec 

QUP  Upstream  inflow  of  water  to  segment  ft3/sec 

RAIN  Rain  precipitat ion  inches 

RNET  Net  rainfall  intensity  in/hr 

S  Specific  gravity  of  sediment 

TC  Critical  bed  tractive  force  lbs/ft2 
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UST  Bed  shear  velocity 

VEL  Average  velocity  of  flow 

VINT  Total  interception  per  unit  area  of  land 

VS  Settling  velocity  in  quiescent  water  for  material 

fraction 

X  Distance  measured  from  upstream  boundary  of  segment 

XS  Distance  travelled  by  a  characteristic  or  shock 

water  wave  during  the  current  time  step 
XSS  Space  increment  used  for  sediment  routing  (Eq.  81) 

ZD  Thickness  of  detached  soil  for  any  one  material 

fraction 
Bed  elevation 


ft/sec 

ft/sec 

inches 

ft/sec 

ft 

ft 

ft 


ZL 


ft 

ft 
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STREAM  CHANNEL  STABILITY.  APPENDIX  l.  SIN8LE  EVENT  NUMERICAL  MO— CTC<U> 
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PROGRAM  SEDLAB 

C  ************** 

c 

DIMENSION  TITLE  <20  >  . SEG < 1 0 0 ) , 1  SEC ( 1 00 ) , QOUT < 300 ) , GTOT ( 30 0 > , 
*SLEN<100> , SLOPE (100) ,CPER ( i 0 0  > , EPER < 100 ) , I AR Y ( 1 00 ,5) , OMESI300 ) , 
♦GMES<300  )  ,ZZ<3> ,C<?1 >  ,OVA<  100)  ,NEED<6)  .CANOUOO)  ,GCOV<  100  >  , 
*HYCND< 100  > ,SORPTY< 100) ,  STORM! 5) 

INTEGER  SEG 

COMMON  /GEN/  UMAX , I TGOM , ITPON,DT ,DTS 

COMMON  /CO V/  CND.GCD ,HLR ,HIR 

COMMON  /RAI/  DR <300) , RNET  <  300  > ,QINF<300) 

COMMON  /INT/  VOR , VOG , SRG , VI N , EVP , VINT , RAIN 
COMMON  /INF/  PERM, SORP 

COMMON  /FLO/  0 <30 0 , S ) , ALP ,EXP , K , ITYPE , SLN , SLP ,CPR ,EPR , 
*KOUT,QUP<300  >  ,(3LAT<300> 

COMMON  /SED/  GS < 30 0 , 5 , 5 ) , ZL < 1 0 0 > , ZD < 1 0 0 , S ) , CDUP <30 0 , S ) , 

♦GDLAT (300,5)  , POR  < S )  , CP  < 5 > , DSO < S ) , PC < S > , VS< 5 > , ERD< 300 ) , 

♦  BEGQOO)  , END! 100  >  ,  ADF  ,  SNU  ,  EMV ,  HOC  ,  SGC  ,  QF. ,  AP  , K 1 , BRE , APS , I TS , X , 
*XS , XSS , NO , BE  T , DX , SUBW , SPGR , ACCLN , GAMA , K  2 , K  S , GMAX , NSO IL 
DATA  ZZ/' I ' , 'M' , 'C '/ ,DEC/' . '/.BLANK/'  '/ 

C 

C  DATA  INPUT 
C 

READ! S , SO  1 )  TITLE 
C  WATERSHED  GEOMETRY 

READ! S, 502)  ARE A, NOV ,HCH ,NCI ,NSTRM 

NTO=NOV+NCH 

NSEG=NTO+NCI 

IC=NOV+l 

READ<  5 ,503 )  <SEG<1) , OVA < I > , SLEN < I ) , SLOPE < I ) , CPER < I ) , EPER < I ) , 
*1=1, NOV) 

READ (5,504)  <  SEG < 1 > ,SL  EN<I ) ,Sl OPE  < I ) , CPER ( I > , EPER < I > , I=IC ,NSEG > 
C  VEGETATIVE  COVER  AND  SOIL  CHARACTERISTICS 
READ<  5  ,  SOS )  VOG  ,SP.G  ,  VOR  ,HLR 
READ!5,S06)  NSOIL.SPGR , AIM, ADF , GMAX, GDX 
READ< 5 ,516)  <DS0  <  I  )  ,  I ---1  , NSOIL  ) 

READ  <5,517)  < PC < I ) , I  =  1 , NSOI L > 

READ <5, 50 7)  <CANO<I) ,CCOV<I) ,HYCND<I), 

*1=1 ,NSEG ) 

C  COMPUTATIONAL  SEQUENCE 

READ<  5,500)  <  1SEC  <  I )  ,  <  IAR  Y  <  I  ,  J  )  ,  J  =  1  ,  S )  ,  I  -- 1  ,  NSEG) 

C  WATER  PROPERTIES  AND  FLOW  ROUTING  PARAMETERS 
READ <5, 509)  TEMP ,GAMA , SNU  ,EXP ,C1 ,C2 
C  OUTPUT  OPTIONS 

READ  <  S , 5 1 0 )  <NEED< I ) , 1  =  1 ,6) 

C 

C  LISTING  OF  INPUT  DATA 
C 

WR ITE <  6 , 60 1 )  TITLF 

IF  <  NEED  < 1 )  EQ  0)  GO  TO  9S 

WRITE  <6, 603)  AREA  ,  NOV  ,  NC.H  ,  NCI  ,  NSEG  ,  NSTRM 

WRITE <6, 604 )  <SEG< I ) ,OVA< I > , SLEN ( I ) , SLOPE < I ) , CPER < I ) ,EPER<I ) , 
*CANO< I ) , GCOV ( I ) ,HYCND  < I ) ,1  =  1 ,NOV> 

WRITE <6 ,605)  <SEG< I ) , SLEN  < I ) , SLOPE! I ) ,CPER<I ) , EPER < I ) , CANO< I ) , 
*GCOV  <  I )  ,  HYCND  <  I  )  ,  I--IC.NSEG) 

WRITE  <6, 606)  VOG  ,  SRC, ,  VOR  ,  HLR  ,  NSOIL  ,  SPGR  ,  AIM ,  ADF  ,  GMAX ,  GDX 

WRITE  <  6 , 6 1 6 ) 

WRITE <6 ,61 7)  <I,D50<I) , PC < I > , I = 1 , NSOIL) 

WRITE <6, 618) 

WRITE <6, 607)  <ISEG< I ) ,SEG<I ) , <IARY(I,J> , J=1 ,5) , 1=1 ,NSEG> 
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URITE(6 ,606)  TEMP , GAMA, SNU, EXP ,Ci,C2 
95  CONTINUE 
C 

C  INVARIANT  INFORMATION 
ACCLN=32 . 2 

SUBW=(SPGR-1 . 0 )*GAMA 
C  POROSITY  AND  SETTLING  VELOCITY 
DO  88  1=1 , NSOIL 

POR ( I ) =1 . -0.245-0 . 0864/(0 . 1*D50 < I ) >**0 . 21 
DMM=D50 ( I > 

CALL  SETVEL<DMM,TEMP ,SPGR ,W) 

VS(I)=W 

D50(I)=D50(I)/304.8 
88  CONTINUE 
C 

TRAIN=0 .0 
TVINT=0 . 0 
C 

C  COMPUTATION  FOR  EACH  STORM 
DO  100  N=1 , NSTRM 
C 

C  STORM  DATA  INPUT 

READ! S, SI 1 )  STORM, DTM , ITMAX , ITCOM , EVP ,VIN 
READ  <  S , S12 )  ( SORP  TY  <I),I  =  1,NSEG) 

READ( S ,513)  <DR(I> ,1=1 , ITMAX) 

READ  <  S , S14 )  <QMES(I ) ,1  =  1, ITCOM) 

READ(S , 515 )  (GMES(I), 1=1, ITCOM) 

DT=DTM/60 . 0 
DTS=DTM*60 .0 
C 

C  LISTING  OF  STORM  DATA 

IF <NEED< 1 )  . GT . 0  )  URITE(6,6ii>  STORM ,VIN , EVP , DTM , ITMAX , ITCOM 
IF(NEED< 1 ) . GT . 0 )  URITE<6,612>  < SORPTY ( I > , 1=1 , NSEG > 

C 

C  POTENTIAL  EROSION  BY  RAINDROP  IMPACT 
DO  130  IT=1, ITCOM 
IFUT.GT.  ITMAX)  DR(IT)  =  0.0 
ERO ( I T ) =AIM* ( DR ( IT ) /4320  0 . 0)**2. 0 
130  CONTINUE 
C 

C  ROUTING  SEGMENTS  ACCORDING  TO  COMPUTATIONAL  SEQUENCE 
UR ITE ( 6 , 621 ) 

IF (NEED<2) . GT . 0 )  WRITE<6,648> 

DO  101  1=1, NSEG 
C 

K=ISEG< I ) 

IF(K.LE.NOV)  ITYPE=i 

IF<  K . GT  NOV . AND . K . LE . NTO )  ITYPE  =  2 

IF  <  K , GT  NTO )  ITYPE=3 

SLN=SLEN<  K ) 

SLP  =SLOPE  <  K ) 

CPR=CPER ( K ) 

EPR=EPER  <K ) 

CND=CANO  <  K  > 

GCD=GCOV ( K ) 

PERM=HYCND(K ) 

SORP=SORPTY (K ) 

EMV=HLR#< 1  0-GCD) 

IF  <  GCD . EQ . 0 , 0 )  EMV=0 . 0 
HIR=0 . 0 
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HGC=HIR*GCD 

5GC=GCD-HGC 

E»RE=<  1  0  -CND )  *  <  1  0  -CCD  ) 

QVF=OUA(K >/<AREA*43560  0) 

C 

C  SELECT  SPATIAL  INCREMENT  FOR  SEDIMENT  ROUTING 
RDX=SLN/GDX 
NDX=INT (RDX ) 

IF ( (RDX-NDX ) . GT . 0 . 0  )  NDX=NDX+1 
IF  <  NDX .  LT . 2 )  NDX  =  2 
DX=SLN/FLOAT (NDX) 

C  KINEMATIC  FRICTION  PARAMETERS 
ALP=C1+C2*SLP**0  5 
C 

C  INITIALIZE  SURFACE  ELEUATION  AND  DETACHED  SOIL  STORAGE 
DO  107  J=i,100 
Zl.<  J)=0 . 0 
DO  87  L=i, NSOIL 
ZD ( J , L  >  =  0  0 
87  CONTINUE 

107  CONTINUE 

IF ( ITYPE . ER . 2)  GO  TO  108 
C 

C  COMPUTING  NET  RAINFALL  RATE 
CALL  INTCPN 

IF ( ITYPE . EQ  i)  TRAIN =TRAJN+RAIN*OUF 
IF  < ITYPE . EQ . 1)  TVINT-TVINT+UIN7#OUF 
C 

C  COMPUTING  PONDING  TIME  AND  POTENTIAL  INFILTRATION 
CALL  INFLTN 

IF  (  UPON  .  EQ  .  0  )  GO  TO  151 

108  CONTINUE 
C 

C  WATER  AND  SEDIMENT  ROUTING 

C 

C  UPSTREAM  INI  LOU  AND  LATERAL  INFLOU/OUTFLOW 
DO  102  IT=ITPON, ITCOM 
QUP ( I T )  =  0 . 0 
QLAT <  IT )==0  0 
DO  81  L= 1 , NSOIL 
GDUP <  I T  ,  L  >  =  0  . 0 

81  GDLAT ( IT  >L ) ~0 .0 

IF< IARY(K ,3) . EQ . 0)  GO  TO  103 
DO  104  J=-l  ,2 

IF(IARY<K ,2+J)  EQ  0)  GO  TO  103 
J  J  =  I AR Y ( K , J+2) 

QUP (IT) =QUP (IT)+Q< IT, JJ> 

DO  82  L=i , NSOIL 

82  GDUP (IT, L) =  GDUP ( I T , L  ) +GS ( I T  ,  J  J , L  ) 

104  CONTINUE 

103  IF< ITYPE. GT  1)  GO  TO  505 

QLAT  < IT )=QLAT ( IT )  +  ( RNET < IT ) -QINF ( IT ) )/4320  0 . 0 

105  IF ( I ARY(K , 1  )  .  EQ . 0  >  GO  TO  102 
DO  106  J  =  1  ,? 

IF(IARY(K,J).EQ.O)  GO  TO  102 
JJ=IARY(K , J ) 

QLAT  (  I T  )  -  Q  L.  A  T  (  IT  )+Q(  IT  ,JJ> 

DO  83  L~1 , NSOIL 

83  GDLAT  ( I T  ,  1. )  -  GDLAT  <  I T  ,  L  )  ■*  GS  <  I T  ,  J  J  ,  L  > 

106  CONTINUE 


102  CONTINUE 

KOUT  =  I  ARY  <  K ,S> 

C 

CALL  UROUT 
C 

C  LISTING  THE  FINAL  SURFACE  LEVEL  OF  THE  SEGMENT 
IF ( NEED ( 2  > . GT . 0  >  URITE(6,649>  K 

IF  <NE£D<2>  . GT  0 )  URITE(6,650>  ( BEG ( J ) ,END < J > , ZL ( J > , J*1 , NO ) 

C 

101  CONTINUE 
C 

C  RUNOFF  AND  SEDIMENT  DISCHARGE  AT  THE  CATCHMENT  OUTLET 
DO  121  I T  *  1 , 1 T COM 
IF< IT . C£  ITPON)  GO  TO  122 
QOUT  (  I  T  >  =  0  0 
GTOT(IT)*0  0 
GO  TO  121 

122  QOUT (  1 T  )  =Q ( I T ,  1  > 

GTOT  < I T  >  =  0  0 

DO  84  L=l,NSOIL 

84  GTOT  < IT ) -GTOT (IT) «CS( I T , 1 , L  >  *GAMA 
121  CONTINUE 

UR  I TE ( 6 , 622 ) 

C 

C  COMPUTED  YIELDS,  PEAKS  AND  TIME  TO  PEAKS  OF  WATER  AND  SEDIMENT 
TQOUT  =0  0 
QCPK=0  0 
TGOUT  =  0  0 
GCPK=0  0 

DO  123  I T  =  ITPON , I TCOM 

IF (QOUT (IT) .LT  .  QCPK )  GO  TO  124 

QCPK=QOUT ( IT ) 

IQC=IT 

124  TQOUT  =  TQOUT  +  QOUT  <IT)#DTS 

IF (GTOT ( IT ) . LT . GCPK )  GO  TO  125 
GCPK=GTOT( IT) 

IGC=IT 

125  TGOUT  = T GOUT + GTOT ( IT ) *DTS 

123  CONTINUE 

TQPSF= TQOUT /(AREA* 43560 . 0>*12  0 
GO  TO  126 
151  CONTINUE 
C 

C  SMALL  EVENT  GENERATING  NO  RUNOFF 
TQOUT=0 .0 
TQPSF=0 . 0 
QCPK  =  0  .  0 
IQC  =  0 
TGOUT=0 . 0 
GCPK=0 . 0 
IGC  =  0 

126  CONTINUE 

IF (NEED<3) .EQ. 0 )  GO  TO  96 
C 

C  WATER  BUDGET 

TINF=TRAIN-(TVINT+TQPSF> 

C 

C  MEASURED  YIELDS,  PEAKS  AND  TIME  TO  PEAKS  OF  WATER  AND  SEDIMENT 
TQMES=0 . 0 
QMPK=0 . 0 
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TGMES=0 . 0 
GMPK- 0  0 

DO  127  IT  =  1  ,  ITCOM 
IF<QM£S(IT)  LT.QMPK)  GO  TO  128 
QMPK=QM£S< IT) 

IQM=IT 

128  TQMES~TQMES+QMES<  IT)X<DTS 
IF(GMES( IT ) . LT . GMPK )  GO  TO  129 
GMPK=GMES( IT) 

IGM=IT 

129  TGMES=TGMES+GMES(IT)*DTS 
127  CONTINUE 

PERCENTAGE  DIFFERENCE  OF  COMPUTED  VALUES  TO  MEASURED  VALUES 
PCWY  =  <  TQMES-TQOUT) /T QMES* l 0 0 . 0 
PCWP=(QMPK~QCPK) /QMPK* 100  .0 
PCWT =FLOAT (IQM-IQC)/FLOAT(IQM>*100  .  0 
PCSY=<TGM£S-TGOUT)  /TGMF.S*  1  0  0  .  0 
PCSP=  <  GMPK-GCPK ) /GhPK$l 0 D . 0 
PCST  =  FLOAT< IGM-IGC > /FLOAT (IGM>*100  .  0 

OUTPUT  OF  WATER  BUDGET 

WRITE ( 6 , 6S1 )  TRAIN,TVINT,TINF,TQPSF 

OUTPUT  OF  YIELDS,  PEAKS  AND  TIMF  TO  PEAKS  OF  WATER  AND  SEDIMENT 

WRITE (6,652)  TQMES , TQOUT , PCWY , QMPK . QCPK , PCWP , IQM , IQC , PCWT , TGMES , 
♦TGOUT ,PCSY ,GMPK ,GCPK ,PCSP , IGM , IGC , PCST 
CONTINUE 

IF (NEED(4) .EQ . 0 )  GO  TO  115 

OUTPUT  OF  RAINFALL  INTENSITIES,  AND  MEASURED  AND  COMPUTED  HYDROGRAPHS 
WRITE  <  6 , 656 ) 

IP=ITPON-l 

WRITE  (6, 657)  <  J  ,  DR  <  J  )  ,RNF.T  <  J  )  ,  QMES  ( J  )  ,QOUT(J>  ,GMES(J)  ,GTOT(  J)  , 

*J=1 , IP ) 

WRITE <6, 658)  < J , DR ( J ) , RNET ( J ) , QINF ( J ) , QMES ( J ) ,QOUT ( J ) , GMES ( J ) , 
*GTOT< J) , J=ITPON, ITCOM) 

115  CONTINUE 
C 

C  OUTPUT  OF  HEITOGRAPH ,  MEASURED  AND  COMPUTED  HYDROGRAPH 
IF ( NEED  <  5  >  .EQ.0)  GO  TO  97 
UR  I TE ( 6 , 66  0  ) 

QMAX=  QMPK 

IF(QCPK .GT .QMPK)  QMAX=QCPK 
DO  141  J=1 ,71 

141  G(J)=BLANK 

DO  142  J" 1 , ITCOM  ,  2 
G( 1  )  =  DEC 
IR=DR ( J ) *7 . 0 
IRR=71  -IR 
DO  143  J1 =IRR , 7 1 

143  G(J1)=ZZ<1> 

IM=QMES< J)/QMAX*70 . 0 
G  < IM  >  =ZZ  <  2 ) 

IC=QOUT< J)/QMAX*70 . 0 
G < IC ) -ZZ ( 3  > 

WR ITE  <  6 , 66 1 ) J , G 
DO  144  Jl-1,71 

144  G  <  Ji ) =BLANK 

142  CONTINUE 


97  CONTINUE 
C 

C  OUTPUT  OF  HEITOGRAPH,  MEASURED  AND  COMPUTED  SEDIMENTGRAPH 
IF <NEED(6 ) . EQ . 0 )  GO  TO  98 
WRITE<6,662> 

GMAX=GMPK 

IF  (GCPK  .  GT  .  GMPX  )  GMAX«=GCPK 
DO  161  J=1 ,71 

161  G< J)=BLANK 

DO  162  J=1,ITC0M,2 
G( 1  )*DEC 
IR*DR<J>*7. 0 
IRR-  71-IR 
DO  163  Ji=IRR ,71 

163  G<J1)=ZZ<1) 

IM=GMES<J)/GMAX*70 . 0 
G ( IM) =ZZ  <  2 ) 

IC=GTOT< J)/GMAX*70 . 0 
G(IC)=ZZ(3> 

WRITE ( 6 , 663  >J , G 
DO  164  Jl=l ,71 

164  G< Ji)=BLANK 

162  CONTINUE 

98  CONTINUE 
WRITE ( 6 ,65? ) 

100  CONTINUE 
C 

SOI  FORMAT  <20A4 ) 

SU2  FORMAT'FiO .4,414) 

503  FORMAT (I4,F12.3,4F10.4) 

504  FORMAT <I4,4F10.4) 

505  FORMAT ( 5F1 0 . 4 ) 

506  FORMAT ( I4,F10 ,4,2F10  8,2F10 .4) 

516  FORMAT ( 1 0F7 . 3  > 

517  FORMAT (10F7. 3) 

507  FORMAT (3F10 . 4 ) 

508  FORMAT  <614 ) 

509  FORMAT ( 2F1 0 .4,F10.8,3F10.4> 

510  FORMAT (614) 

511  FORMAT (SA4 ,F 1 0 . 4 , 2 1 4 , 2F 10.4) 

512  FORMAT  <  8F7 . 4 ) 

513  FORMAT(12F6.3) 

514  FORMAT ( 12F6 . 2 ) 

515  FORMAT  < 12F6 . 2 ) 

601  FORMAT (1H1///X,20A4///) 

603  FORMAT(30X, 'CATCHMENT  GEOMETR Y ' //4 IX , '  DRAINAGE  AREA  =  ', 

♦F7.1,'  ACRES  '/38X, 'OVERLAND  SEGMENTS  =  ', IS/10X ,' CHANNEL  ', 
♦'SEGMENTS  WITH  NEGLIGIBLE  INFILTRATION  =  ', I5/8X CHANNEL  ', 
♦'SEGMENTS  WITH  CONSIDERABLE  INFILTRATION  =  ', 15/41 X ,' TOTAL  ', 
♦'SEGMENTS  =  ', IS//30X , 'NUMBER  OF  COMPUTED  STORMS  =  ',15//// 

♦24X, 'PHYSICAL  CHARACTERISTICS  OF  SEGMENTS'//  X,' - 

♦ - - - ' 

♦/36X, 'RELATION  BETWEEN  CANDPY  GROUND  SATURATED ' /X ,' SEGMENT  OVE 
♦RLAND  LENGTH  SLOPE  WETTED  PERIMETER  COVER  COVER  HYDRAULIC 
♦'/X, 'NUMBER  AREA' ,21X, 'AND  FLOW  AREA  DENS.  DENSITY  CONDU 

♦CT.'/X,'4  TYPE  <SQ  FT)  <FT>  < FT/FT )  COEFF .  EXP.',20X, 

♦  '(  INCH/HR  >'/X,  ' - ' 

♦  ,' - '//) 

604  FORMAT < X , 13 , '  OV ' , X , F 1 0 . 1 , X , F7 . 1 ,  2X ,F7 . 5 , 2X ,F5 . 1 , 3X , F5 . 3 , 4X , F5 . 3 , 
♦3X,F5.3,4X,F5 .3) 


605  FORMAT <X, 13, '  CH' , 13X ,F7  1 ,X ,F7 . S ,2X ,FS . 1 ,3X ,FS . 3 , 4X ,FS . 3 ,3X , 
tF5.3^X,FS  3) 

606  FORMAT < ////27X , 'GROUND  COVER  CHARACTER  1ST  ICS ' //6X , 

♦  ' INTERCEPTION  STORAGE  CAPACITY  FOR  GROUND  COVER  ■=  ', 

♦F9  3,'  INCHES' ,/6X, 'RATIO  OF  EVAPORATING  SURFACE', 

♦  '  TO  PROJECTED  AREA  =  ' , F9 . 3/7X , ' R AT  10  BETWEEN  INTERCEPTION'. 

♦  '  STORAGE  CAPACITIESV19X  ,  'OF  CANOPY  COVER  AND  GROUND  COVER', 

♦  '  -  ',F9.3/10X, 'AVERAGE  HEIGHT  OF  GROUND  COVER  IN  CHANNELS', 

♦  '  =  '  ,F9.3,  '  FEET' 

♦  Z///18X, 'SOIL  CHARACTERISTICS  AND  DETACHMENT  PARAMETERS'// 
♦27X, 'NUMBER  OF  SEDIMENT  SIZE  FRACTIONS  =  ',15/ 

♦32X, 'SPECIFIC  GRAVITY  OF  SEDIMENT  -  ' ,FB.2/liX, 

* 'COEFFICIENT  OF  SOIL  DETACHMENT  BY  RAINDROP  IMPACT  =  ', 

♦E10  3/17X, 'COEFFICIENT  OF  SOIL  EROSION  BY  SURFACE  FLOW  = 

♦E10 . 3/22X, 'MAXIMUN  PENETRATION  DEPTH  UF  RAINDROPS  =  ',F9.3, 

♦  '  FEET '/20X , 'SPACE  INCREMENT  USED  IN  SEDIMENT  ROUTING', 

♦  '  =  ' ,F8.2,'FEET'//> 

616  FORMAT  < 1SX , ' - 

♦/16X, 'SERIAL  NO.  OF  MEDIUM  SIZE  OF  PERCENTAGE ' /I 6X , 

♦  'SED  FRACTION  SEDIMENT  FRACTION  OF  SEDIMENT ' /37X , 

♦ '  <MM)'  ,13X,  'FRACTI0N'/15X,  ' - 

♦  ' -  - . /r> 

617  FORMAT (19X,I4,l3X,F7.3,12X,FS.i> 

618  FORMAT  <////SX,  ' -  - 

♦  ,' - '/9X,  'ISEG'  ,30X,  '  IARYV27X, 'THIS  ARRAY  ' 

♦'CONTAINS  THE  STORAGE  LOCATI ONS ' /SX , ' COMPUTAT IONAL ' , 9X , ' OF  ', 
♦'COMPUTED  INFLOWS  AND  OUTFLOWS  IN  THE ' /8X ,' SEQUENCE 1 IX ,' WATER ' 

♦  '<(■>)  AND  SEDIMENT  <  GS>  DISCHARGE  MATRICES'/ 

♦  22X  ,  ' - '/ 

♦  11X, '♦', 1  OX , 'SEGMENT  t  LATERAL  INFLOWS  UPSTREAM  INFLOWS  ', 

♦  'OUTFLOU'/SX,' -  -  -  ', 

♦  ' -  - '//) 

607  FORMAT <1  OX, 13 , 1 1 X , 13 , 9X , 1 1 , 6X , 1 1 , 9X , 1 1 , 8X , 1 1 , BX , 1 1 > 

608  F0RMAT<////18X, 'WATER  PROPERTIES  AND  KINEMATIC  FRICTION 
♦ ' P AR AMETER S ' //3 1 X , 'WATER  TEMPERATURE  =' 

♦  ,F8 . 2  ,  '  I)EG  CENT  '/24X, 'SPECIFIC  WEIGHT  OF  WATER  *' 

♦,F8.2,'  LBS/CU.KT  '/20X, 'KINEMATIC  VISCOSITY  OF  WATER  =  ', 

♦  HX  ,  El  0 . 3 , '  SQ.FT. /SEC. '/13X, 

♦'EXPONENT  IN  KINEMATIC  WAVE  EQUATION  =  ' ,F7  2/19X, 

♦'KINEMATIC  FRICTION  PARAMETERS  ,  C1=',F6.2,',  C2  =',F6.2//> 

611  FORMAT ( //36X, 'STORM  DATA ' //37X , ' STORM  DATE  =  ',5A4/16X, 

♦  'ANTECEDENT  INTF  RCEPTION  STORAGE  =  ',F10.4,'  INCHES'/26X, 

♦'MEAN  EVAPORATION  RATE  =  ',F10.4,'  INCHES  PER  HOUR'// 

♦  2SX , ' T I  ME  STEP  DURING  STORM  *',F8.1,'  MINUTES'/ 

♦8X , ' NUMBER  OF  TIME  STEPS  IN  RAINFALL  PERIOD  =  ', 

♦I5/6X, 'NUMBER  OF  TIME  STEPS  IN  SIMULATION  PERIOD  =  ', 

♦I5///21X, 'SORPTIVITY  FOR  EACH  CATCHMENT  SEGMENT', 

♦/28X, ' (SQUARE  INCHES  PER  HOUR)'/) 

612  FORMAT ( X , 1  OF  8 . 4  > 

621  FORMAT ( ////20X , ' ***###***#*****##* #*#***#*********####*##** ' / 

♦20X,'*  START  OF  WATER  AND  SEDIMENT  ROUTING  *' 

♦  /20X,  ' *****#********#*###*##**#*#**) ***##****#**#*') 

648  FORMAT<////?X,  ' - '  , 

%  ' - ' /9X  'DISTANCE' 

♦'  FROM  U/S  END  DISTANCE  FROM  U/S  END  CHANGE  OF  BED  ELEV.'/ 
♦8X,'0F  SEGMENT  TO  START  OF  OF  SEGMENT  TO  END  OF  ' , 

♦'FOR  THE  SPACE'/iOX, 'SPACE  INCREMENT  ( FT )', 2X , ' SPACE 

♦  'INCREMENT  ( FT ) ' , SX , ' I NCREMFNT < FT ) '  /9X  ,  ' - 

♦  ' -  - //> 

649  F0RMAT<36X, 'FOR  SEGMENT  ',14) 


650  FORMAT  (1SX.F7. 1  ,  17X ,F7  1  ,  1SX ,F6 . 3 > 

622  FORMAT  <////15X , ' a******************************************* ' , 

♦'*»*#***'/15X,  '*  END  OF  WATER  AND  SEDIMENT  ROUTING', 

♦  '  *'/15X, '*******************»#********«' , 

S'**********#**##*****#'//) 

651  F0RMAT(//35X, 'WATER  BUDGET ' //32X , 'PRECIPITATION  =  ',F8.2, 

♦  '  INCHES'/13X, 'LOSS  DUE  TO  INTERCEPTION  STORAGE  =  ',F8.2, 

♦  '  INCHESV21X, 'LOSS  DUE  TO  INFILTRATION  -  ',F8.2,'  INCHES'/ 

♦34X, 'WATER  YIELD  =  ',F8.2,'  INCHES'///) 

652  FORMAT <//33X,  'SUMMARY  OF  RESULTS '/X,' - ', 

♦  ' - '  , 

♦  ' - '/8X, 'ITEM  ' , 1SX, 'MEASURED' ,7X, 'COMPUTED' ,9X, 'UNITS' ,8X, 

♦  'PER  CENT  '  /73X  ,  'DIFFER  .  '/X,  ' - '  ,2X,  ' - '  , 

♦  ' - '  ,2X, ' -  -  - '// 

♦  6X, 'WATER  YIELD' ,7X,F12. 1 ,3X,F12 . 1 ,6X, 'CUBIC  FEET ' ,6X ,F6  1 , /X , 

♦ 'WATER  DISCHARGE  PEAK'  , 7X , F8  1 , 7X , F8 . 1 , 4X, ' CUBIC  FEET/SEC', 
♦4X,F6.i/2X, 'TIME  TO  WATER  P  F.AK  '  ,  9X  ,  1 5 , 1  OX ,  IS  ,  9X ,  '  T IME  STEPS', 
♦5X,F6. 1/4X, 'SEDIMENT  YIELD ' , 6X , F 1 2 .  1 , 3X , F 12 . 1 , 8X , ' POUNDS ' ,8X, 
♦F6.1/2X,'SEDIM.  DISCHARGE  PEAK ', 5X , FB . 1 , 7X , F8 . 1 , 5X ,' POUNDS  ', 
♦'PER  SEC',  F9 . 1 /X , 'TIME  TO  SEDIMENT  PEAK ' ,7X , 15, 10X , IS , 

♦9X  ,  'TIME  STEPS' ,SX,F6 . l//> 

656  FORMAT  ( ///X,  ' - ', 

♦ ' - '/X,  'TIME'  ,3X  , '  PRECIP  ' 

♦,SX, 'NET' ,6X, 'POTENTIAL' ,9X, 'RUNOFF' ,8X, 'SEDIMENT  DISCHARGE ' /X , 

♦  'STEP' ,11X, 'RAINFALL  INFILTR  AT  ION '  ,  2X  ,  ' - ' 

♦  ,2X ,  ' - '/40X,  'MEASURED  COMPUTED  MEASURED  ', 

♦'COMPUTED' /3X, '#' , SX , ' IN/HR ' , 3X , ' IN . /HR . ' , SX , 'IN. /HR  '  ,7X, 

♦'CFS' ,7X, 'CFS' ,5X, 'LBS/SEC' , 3X , ' LBS/SEC ' /X , ' -  -  ', 

♦  -  -  -  -  -  - ' 

♦  /) 

657  FORMAT <X, I A,3X,F6  3 ,3X ,F6 . 3 ,8X , '-' ,BX,F7. 1 ,3X,F7. 1 , 

♦3X,F7  1,3X,F7.1) 

658  FORMAT <X,I4,3X,F6.3, 3X ,F6.3,6X,F6.3 ,5X , F7 . 1 , 3X ,F7 . 1 , 

♦3X,F7 . 1 ,3X,F7 . 1 ) 

659  F0RMAT<//8< ' ********** ' ) ) 

660  F0RMAT<///4X,7( '♦****♦****' ), '*'/4X, '*  PLOTS  OF  HYETOGRAPH<  I )  ■,  ' 

♦  '  COMPUTED < C )  AND  MEASURED <M>  HYDROGRAPHS  *'  /4X, 


♦7( '****»**Ht**' ) , '*'//14X , ' (Q/QMAX)*i00  - >',19X,'< -  PRECIP', 

♦' (INCHES/HOUR) '/2X, 'TIME'/2X, 'STEP ' ,37X, '5  4  3 

♦'  1  0'/3X,'#',4X,'0  10  20  30  40 

♦  '  60  70  80  90  10  0'  /2X ,  ' -  I',10<' - 1')) 

661  FORMAT (2X,I4,2X,71A1 ) 

662  F0RMAK///4X  ,7<  '*****#****'),'#' /4X,  '*  PLOTS  OF  HYETOGR APH<  I )  ;  ' 

♦  '  COMPUTED < C )  AND  MEASURED ( M >  SEDIGRAPHS  *'  /4X, 

♦  7<  '*****#**##'),  ')K'//14X,  '  <G/GMAX)!K100 - >',19X,'< - PRECIP', 

♦ ' ( INCHES /HOUR )'/2X, 'TIME'/2X, 'STEP', 37X , ' 5  4  3  2', 

♦  '  1  0'/3X, '♦' ,4X, '0  10  20  30  40  50' 

♦  '  60  70  80  90  100'  /2X, ' -  I',10<' - 1')> 


663  FORMAT  <2X,I4,2X,71A1) 

END  FILE  6 

STOP 

END 

SUBROUTINE  INTCPN 

C  THIS  SUBROUTINE  COMPUTES  INTERCEPTION  LOSSES  AND  NET  RAINFALL 
C  INTENSITY 

COMMON  /GEN/  ITMAX , ITCOM , ITPON ,DT ,DTS 

COMMON  /COV/  CND , GCD , HLR , HIR 

COMMON  /RAI/  DR(300) ,RNET<300) ,QINF<300> 

COMMON  /INT/  VOR,UOG,SRG,VIN, EVP, UINT, RAIN 
C 


on  ooooo  oro  n  non 


C  INITIALIZING. CUMMULATIVE  VALUES  AND  STORAGE  CAPACITIES 
TR1=0 . 0 
TR2=0 . 0 
TR3=0 . 0 
TR4=0 .0 
TRC=0 . 0 
TRG=0 . 0 
VOC=VOR*VOG 
SRC=VOR*SRG 
CCAP  =  <  I  .  0-VIN)*VOC 
GCAP=(i . 0-VIN)*VOG 

COMPUTING  NET  RAINFALL  FOR  EACH  TIME  INTERVAL 
DO  201  IT=1  ,  ITCOM 

RATE  OF  INTERCEPTION  ON  CANOPY  AND  THROUGHFALL  FROM  CANOPY 
TR1=TR1+DR(IT)*DT 
IF  <  TRi . LE . CCAP )  DRC  =  DR  ( IT  > 

IF<  TR1  .  GT  .  CCAP  >  DRC=--EVP*SRC 

IF (TRi  GT . CCAP  AND  TRC.LT. CCAP )  DRC= ( CCAP -TRC ) /DT 
TRC  =  TRC-»DRC*DT 
RTHO=DR< IT>-CND#DRC 

RATE  OF  INTERCEPTION  ON  GROUND  COVER  AND  NET  RAINFALL  RATE 
TR2=TR2+RTH0*DT 
IF  <  TR2 . LE . GCAP  >  DRG  =  RTHD 
IF<TR2.GT.GCAP>  DRG=EVP*SRG 

IF(TR2  .GT  CCAP . AND . TRG . LT . GCAP  >  DR G~ (GCAP -TRG ) /DT 
TRG--TRG+DRG*DT 
RNET ( IT ) =RTHO-GOD*DRG 
TR3=TR3+RNET  < IT ) *D1 
01  CONTINUE 

TOTAL  RAINFALL  AND  TOTAL  INTERCEPTION  STORAGE 
RAIN-TR 1 
VINT=TRi-TR3 
RETURN 
END 

SUBROUTINE  INFLTN 

THIS  SUBROUTINE  COMPUTES  THE  PONDING  TIME  AND  THE  DECAY  OF 
POTENTIAL  INFILTRATION  FROM  THE  TIME  OF  PONDING 

COMMON  /GEN/  I TMAX , I TCOM , I TPON , DT , DTS 
COMMON  /R AI /  DR(300> ,RNET<300) ,QINF(300> 

COMMON  /INF/  PERM, SORP 

PONDING  TIME 
SL=0 . 0 
ITPON=0 

DO  251  IT  =  1  , ITCOM 
IF  < ITPON . GT . 0  >  GO  TO  252 
SL=SL+RNET ( IT)#DT 
IF(RNET(IT)  LE.PERM)  GO  TO  251 
SR -SORP/ (RNET(IT)-PERM) 

IF(SL.LT.SR)  GO  TO  251 
ITPON-IT 
RP=RNET( IT) 

FP“SL 
AIN=RP 
C=RP-PERM 
C1=C/PERM 
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C2=C/RP 
GO  TO  253 
C 

C  POTENTIAL  INFILTRATION,  THE  DECAY  EQUATION  IS  SOLVED  BY  NEWTON'S 
C  METHOD 

252  SATI=PERM#FLOAT  <IT-ITPON)*DT 
ITR  =  0 

2S A  IF< AIN . LE . PERM)  GO  TO  253 

C3=C2*AIN/(AIN-PERM) 

FF=SATI-FP*<  <RP-AIN>/< AI N-PERM > -Cl *ALOG < C3) ) 

DFF=FP*C*<1  0-<AIN-PERM)/AIN)/CAIN-PERM)**2  0 

AINN=AIN-FF/DFF 

TOL=ABS<  <AINN-AIN)/AINN> 

AIN=AINN 

ITR=ITR+1 

IFOTR.GT.  100)  GO  TO  25S 
IF(TOL . GT . 0 . Oi )  GO  TO  254 

253  IF< AIN . LE . PERM)  AIN=PERM 
QINF< IT)=AIN 

2S1  CONTINUE 
GO  TO  256 

255  WRITE  <  6 , 261 ) 

STOP 

256  RETURN 

261  FORMAT <20X, '**********  ITERATION  EXCEEDS  100  >  INFILTRATION', 
♦  '  EQUATION  **********' ) 

END 

C 

SUBROUTINE  UROIIT 
C 

C  THIS  SUBROUTINE  ROUTES  WATER  FROM  OVERLAND  AND  CHANNEL  UNITS 
C  THE  ROUTING  PROCEDURE  IS  BASED  ON  THE  CHARACTERISTIC  SOLUTION 
C  OF  THE  KINEMATIC  WAVE  APPROXIMATION 

COMMON  /GEN/  ITMAX , IT COM, I  TP ON , DT , DTS 

COMMON  /COV/  CND , GCD , MLR , HIR 

COMMON  /RAI/  DR ( 30 0 ) , RNET ( 30 0 ) , QINF ( 30 0 ) 

COMMON  /FLO/  Q < 30 0 , 5 ) , Al  P , EXP , K , I  TYPE , SLN , SLP , CPR  ,  EPR  , 

♦KOUT , QUP  <300),QLAT<300) 

COMMON  /SED/  GS( 30 0 , 5 , 5)  ,ZL ( 1 0 0  )  , ZD<  1 0  0  ,5 )  ,GI)UP  ( 30  0  , 5  )  , 

♦GDLAT (300,5) ,POR (5) ,CP<5) ,  D50 (5) ,PCC5) , VS (5) ,ERO(300)  , 

♦BEG (100 ) ,END( 100 ) , ADF , SNU , EMV , HGC , SGC , QE , AP ,K 1 , BRE , APS , ITS , X , 
♦XS,XSS,N0,BET,DX,SUBW,SPCR,ACCLN,GAMA,K2,KS,GMAX,NS0IL 
C 

EXP  1 =£XP  +  i . 0 
BET=1 . 0/EXP 
EXM1=EXP-1 . 0 
TERM=EXP*ALP*DTS 
C 

K l=ITPON 
K2"0 

301  IF (K 1 . GT . ITCOM)  GO  TO  391 
KS=K  1 
ITS=K 1 
NO=0 

QB-QUP  <  K 1 ) 

AB- ( QB/ALP  >##BET 

QP  =  QB 

AP=AB 

APS=AP 

DO  381  L=1 , NSOIL 


CP<L>=0  O' 

381  IF ( QP  GT.O.O)  CP <L ) ~GDUP ( K 1 , L > /QP 
QL=QLAT (K1 ) 

IF  ( ITYPE  .  NE  .  3)  GO  TO  302 
APF=AP+QL*DTS 

IF  <  APF  .  GT  .0.0)  PERIM=CPR*APF**EPR 
IF( APF . LE . 0 . 0 )  PERIM=0 . 0 
QL=QL-PERIM*QINF <K1 1/43200 . 0 
302  IF<QL . LE . 0 . 0 . AND . QB .EQ . 0 . 0)  GO  TO  303 
X=0  .0 
XSS-0 . 0 

I F  <  K2 . GT . 0  >  GO  TO  304 
AP=0 . 0 
QP=0 . 0 

DO  382  L=1 ,NSOIL 

382  CP  <  L )  =  0 . 0 
GO  TO  30S 

304  IF (QB . EQ . 0 . 0 )  GO  TO  305 
C  TESTING  FOR  SHOCK  FORMATION 
QA=0 . 0 

IF(Ki.GT.l)  QA=QUP (K 1-1 > 

AA=  <  QA/ALP ) **BET 
IF  < AA . GE . A B)  GO  TO  30S 
SMOCK -AB-AA 
ADD=QL*DTS*0  5 

IF  <K1 . GT . 1 )  ADD-ADD+QLAT(Ki-i )*DTS*0 .5 
IF(SHOCK . LE . ADD)  GO  TO  305 
C 

C  SOLUTION  WITH  SHOCK 
311  IF<K1 . GT . ITCOM)  GO  TO  312 
QL- QLAT  (K  1  ) 

IF  < ITYPE . NT  3)  CO  TO  313 
APF=AP+QL*D1S 

IF < APF  GT  0.0)  PERIM=CPR#APF**EPR 
IF <  APF  .  LE  .  0  .  0  >  PERI M==0  .  0 
QL=QL~PER IM*QINF(Ki)/43200 . 0 

313  ABF=AB+QL*DTS 

IF (ABF .IE  0.0)  ARF  =  0  .  0 

IF ( ABF  EQ . 0  0 )  GO  TO  305 

AAF=AA  +  QL*I)TS 

IF  <  AAF . LE . 0 . 0 )  A AF  =  0  .  0 

QBF=ALP*ABF**EXP 

QAF=ALP*AAF**EXP 

IF ( QL . EQ . 0  0)  GO  TO  314 

DEN=EXP1*<  AB-AA )*QL 

PROD=ALP/I)EN 

XS=PROD*(ABF**EXPl-AAF**EXPl-AB**EXPi+AA**EXPl ) 
GO  TO  322 

314  XS=<QB-QA)*DTS/<AB-AA) 

322  AB=ABF 

AA=AAF 

QB-QBF 

QA-QAF 

X=X+XS 

XSS=XSS+XS 

IF  ( X .  GE  .  SLN  >  GO  TO  32.3 
C 

C  SEDIMENT  ROUTING 

QE  =  ALP*C  (  (QB/ALP  )**BET+( QA/ALP  >#*BET>/2 .  0>**F.XP 
IF(XSS.GT.DX)  CALL  SROUT 
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QP  =  QB 

AP=<AA+A&>/2 . 0 
K 1=K 1+ 1 
GO  TO  311 

323  QB  =  QF<-QL*<X-SLN) 

QA=QA-QL*(X-SLN> 

QC-ALP  *  < <  <  QB/ALP ) **BET+  < QA/ALP  > **BET > /2 . 0  >**EXP 
IF(K1 . EQ  KC)  QC=(GC+Q(KC ,KOUT > )/2 . 0 
Q ( K 1 , KOUT ) =QC 
C 

QE=Q  <  K 1 ,  KOUT ) 

CALL  SROUT 
DO  390  L= 1 , NSOIL 
390  GS(K1 ,KOUT,L)=CP<L)*QE 
GO  TO  324 
C 

C  FLOW  CEASES 
C 

303  IF(K2  EQ.O)  GO  TO  325 
IF <Ki .  GT . KC )  GO  TO  326 
GO  TO  327 

326  K2=  0 

325  Q (K 1 , KOUT ) ~0 . 0 

DO  392  L=1 , NSOIL 

392  GS  ( K 1 ,  K  CUT ,L)=0 . 0 

327  K1=KS+1 
GO  TO  301 

C 

C  SOLUTION  WITHOUT  SHOCK 
30S  IF(K1 . GT . ITCOM)  GO  TO  312 
QL  =  QLAT i K 1 ) 

IFdTYPE.NF.  3)  CO  TO  331 
APF=AP+QL*DTS 

IF  (APE'  GT  0  0)  PERIMd:PR*ArE**EPR 
IF ( APF  t E  0 . 0)  PER  I H- 0  0 
QL=QL-PERIM*QINF<K1  )/43200 . 0 

331  AC  =  AP  +fj)L  #DTS 

I F  <  K  1  FQ.KS  AND  K2  GT  0>  AC=AP+QL*DTS/2 . 0 

IF< AC . I C  0.0)  AC -0  0 

QC=ALP#AC**EXP 

IFIQL.EQ  0.0)  CO  TO  332 

XS=  <  QC-QP ) /QL 

GO  TO  333 

332  XS=TERM*AP**EXM1 

333  X=X+XS 
XSS-XSS+XS 

IF (X  CF . 5LN  AND  K2  EQ.O)  GO  TO  351 
IF ( X . GE . SLN )  GO  TO  334 
IF ( K2  CT  0)  GO  TO  360 
C  INITIAL  RISING  PART  OF  HYDROGRAPH 
Q<K 1 , KOU T  >  =  (QC  +  QP ) /2 . 0 
QE=Q ( K 1  .KOUT) 

CAl  L  SROUT 
DO  393  t  =1 , N 5 O  T  L 

393  G  S ( K 1  ,KOU1  ,L)=CP<L)*QE 
GO  TO  362 

360  CONTINUE 
C 

C  SEDfPFNT  ROOTING 
QE  =  QC 


362 


IFfXSS.GT. DX>  CALL  SROUT 
QP=QC 
AP  =AC 
Kl=Ki+i 
GO  TO  305 
334  Q ( K 1 ,  KOUT  >  =QC-QL# < X-SLN ) 

C 

Q£=Q<K1 ,KOUT> 

CALL  SROUT 
DO  394  L  =  1 , NSOIL 

394  GS<Ki,KOUT,L>=CP(L)*QE 
C 

C  INTERPOLATION  OF  FLOW  AND  SEDIMENT  DISCHARGES  AT  THE  SEGMENT  OUTLET 
324  K3=K i-KC 

IF < K 3 . LE . 1)  GO  TO  3Si 
QD=Q(K1 ,KOUT)~Q(KC,K  OUT ) 

QAD=QD/ r  LClAT  (K3> 

K3M=K3-i 
DO  352  J  =  1 ,K3M 

Q ( KC+ J , K  OUT )  =  Q  < KC ,KOUT  )+QAD*FLOAT<J) 

352  CONTINUE 

DO  395  L~  1 , NSOIL 

GSD=GS < K  1  ,KOUT,L)-GS(KC,KOUT ,L  > 

GSAD=GSD/F1.0AT<K3> 

DO  396  J=1 ,K3M 

396  GS(KC+J ,KGUT , L ) =GS (KC , KOUT , L ) +GSAD4FL0AT ( J > 

395  CONTINUE 

351  KC  =  K 1 

1F(K2.EQ  0)  KC~K 1-1 
IF <KC . LT . 1 )  K C= 1 
K1>=KS+1 

IF ( K2 . EQ . 0 )  K  i~KS 
K2=KS 
GO  TO  301 
C 

C  EXTRAPOLATION  OF  OUTFLOWS  AT  END  OF  LAST  TIME  STEP 
312  IFtKC.EQ  TTCOM)  GO  TO  391 
K3~ ITCOM-K  C 

QD=Q( KC , KOUT  > -Q < KC- 1 ;KOUT> 

QAD-'QD/FLOAT  <  K  3  ) 

DO  353  J  =  1  ,K3 

Q(KC+ J , KOUT ) =Q  <  KC ,  KOUT ) +QAD#FLOAT ( J ) 

353  CONTINUE 

DO  397  L=l, NSOIL 

GSD=GS  <  K  C , K  OUT ,L)-GS(KC-i ,  KOUT ,  L ) 

GSAD=GSD/F'LOAT  <  K3 ) 

DO  390  J  =  1  , K  3 

398  GS < KC+J, KOUT, L)=GS<KC,K OUT, L>+GSAD*FLOAT< J) 

397  CONTINUE 
391  RETURN 

END 

C 

SUBROUTINE  SROUT 

C  THIS  SUBROUTINE  ROUTES  ALL  SEDIMENT  FRACTIONS 
DIMENSION  ADP<5) ,GLAT<5> ,CONC(S> 

COMMON  /GEN/  I TMAX , I TCOM , I TPON  ,  DT , DTS 

COMMON  /COV/  CND  ,  GCI) ,  HLR  ,  HIR 

COMMON  /R AT /  DR < 30 0  >  , RNET < 30 0 > , QINF < 30 0 > 

COMMON  /FLO/  Q < 30 0 , 5 ) , ALP , EXP , K , I T YPE , SLN , SLP , CPR ,EPR , 

*KOUT,QUP <300 > , QLAT  <300  ) 
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COMMON  /SED/  GS(300  ,5 ,5)  ,ZL<  100  )  ,ZD<  100 ,5)  ,GDUP  <  300 ,5  )  , 
*GDLAT<300 ,S> ,POR(S) , CP < S > , D50 ( 5 ) , PC< 5 ) , VS < S ) ,ERO<300> , 
♦BEGtlOO) , END (100) , ADF , SNU , EMV ,HGC , SGC , QE , AP ,K 1 , PRE , APS , ITS ,X , 
$XS , XSS , NO , BET , DX , SUBW , SPGR , ACCLN , GAMA ,  K2 , K  S , GNAX , NSOIL 
C 

C  SOIL  SURFACE  ELEVATION  AND  DETACHED  SOIL  VOLUME  AT  THE  START  OF 
C  THE  CURRENT  TIME  STEP 

IF<K2.EQ.O  OR  K2.EQ.KS)  GO  TO  401 
DO  402  J=l,100 

IF ( X . GE . BEG<  J ) . AND . X . LE . END( J ) >  GO  TO  403 
GO  TO  402 
403  Z=ZL  <  J ) 

DO  490  L=l, NSOIL 

490  ADP ( L ) =  ZD  <  J , L ) 

GO  TO  480 

402  CONTINUE 
GO  TO  480 

401  IF ( K 1 . EQ . KS  >  Z  =  0 . 0 

IF  <K1 . GT . KS)  Z=ZL ( NO ) 

DO  496  L= 1 , NSOIL 
IF(Kl.EQ.KS)  ADP  <  L  )  =  0 . 0 
496  IF ( K 1 . GT . KS )  ADP ( L > =ZD( NO ,L ) 

480  NO=NO+l 

DEG  <  NO ) =X-XSS 
END ( NO ) =X 

C  HYDRAULIC  PARAMETERS 

IF(QE.LT.l.OE-S)  GO  TO  404 
AE=  <  QE/ALP ) *#BET 
VE1.=QE/AE 
UEP=CPR*AE**£PR 
HYR=AE/WEP 
RN=QE/<SNU*WEP  ) 

C  AVERAGE  LATERAL  INFLOW  AND  POTENTIAL  EROSION  RATES  DURING  THE 
C  CURRENT  TIME  STEP 
K4=Ki-ITS+l 
DO  491  L--1, NSOIL 
TGLAT=0 . 0 
DO  4  OS  J  =  I TS  ,  K 1 
4  OS  TGLAT  =  TGLAT  +  GDLAT< J ,L> 

491  GLAT  <L ) =TGLAT/FLOAT  <  K4 ) 

TFRO=0  0 

DO  406  J=ITS,K1 
406  TERO”TERO+ERO<  J ) 

AERO-TERO/FLOAT  <K4) 

DTSS=DTS*FLOAT  <  K  4 ) 

IF(X  GT.SLN)  DTSS=DTS/(1 . + ( X-SLN > / < SLN-X+XS > > 

DIB=BRE*AERO 

C  COMPUTE  EFFECTIVE  TRACTIVE  FORCE  IN  VEGETATED  REACHES 
IF< ITYPE  EO . 1  OR . EMV . Eq . 0 . 0)  GO  TO  408 
IF(Z  GT.EMV)  GO  TO  409 
RATIO=l . 0-Z/EMV 
IF (RATIO . GT . 1 . 0 )  RATIO=i . 0 
EHT=HLR*RATIO 
E5C=SGC*RATI0 
EGC=i . 0-HGC-ESC 

IF  <  HYR . GT . EHT )  EGC  =  1 . 0-HGC-ESC*EHT/HYR 
GO  TO  410 
409  EGC= 1 . 0 -HOC 

GO  TO  410 
408  EGC=i . 0-GCD 
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410  TAO=GAMA*HYR*SLP*EGC 
UST  =  SQR  T  C  TACI4ACCLN/GAMA  > 

DO  499  L=i  , NSOIL 

C  SEDIMENT  CARRYING  CAPACITY 
DFT=D50  <L ) 

DMM=DFT*304  8 
W=VS<L> 

IF { ITYPE . EG) . 1 ) CALL  YALIN ( UST , DFT , SNU , SUBW , SPGR , ACCLN , VEL , HYR , 
♦CEE) 

IFIITYPE.GT  i.AND.DMM.GE.O  1)  CALL  Y ANG ( DFT , UST ,SNU , VEL ,SLP , 
tU, SPGR , CEE) 

IF ( ITYPE . GT . 1 . AND . DhM .  LT  .  0  .  1 )  CALL  LAURSN ( DFT , UST , SPGR  , VEL , 
♦GAMA , SUBW , SNU , DFT , HYR , ACCLN , W , CEE ) 

CONC(L)=CEL*PC(L>/100  .  0 

C  POTENTIAL  SEDIMENT  EROSI ON ( +VE > /DEPOSI TIONC-VE > 

GDPL=< AE*CONC<L)-APS*CP  <L> > /DTSS-GLAT (L > 

BT=W*DTSS/HYR 

IF ( GDPL . GT . 0 . 0  >CO  TO  407 

IF<  BT . LT  .  1 . 0)GDPL=BT*GDPL 

C  DEPTH  OK  DETACHED  SOIL  AT  THE  START  OF  THE  CURRENT  TIME  STEP 
407  ADC=ADP ( L  ) 

C  EROSION  BY  RAINDROP  IMPACT 
RATIO=l . O-ADC/GMAX 
IF  (RATIO  .  l.E  .0.0)  RATIO=0  .  0 

IF(RN  LE  900 . 0)  ADC=ADC+PC < L ) / i 00 . 0*DIB#R AT104DTSS 
RAD=ADC-CDPL/WEP*DTSS 
C  EROSION  BY  SURFACE  RUNOFF 

IF (RAD .GE . 0 . 0  >  GO  TO  411 

DGD=~ADF*RAB 

GD-UEP#  <  ADC+DGD ) /DTSS 

AS=APS*CP<L)  +  (GD+GLAT<L>  >*DTSS 

CGNC(L)=AS/AE 

2D (NO ,  L  )  -  0  0 

DZ=-DGD 

GO  TO  413 

C  DETACHED  SOIL  DEPTH  AND  SOIL  SURFACE  ELEVATION  AT  THE  END  OF  THE 
C  CURRENT  TIME  STEP 

411  ZD  <  NO , L ) -R AD 

DZ=  ( R  AD- ADP  CD)  /F.GC 
413  ZL(NO)-Z  +  I)7/POR  (L) 

412  CP(L)=CONC(L) 

499  CONTINUE 

GO  TO  497 
C  RUNOFF  CEASED 
404  AE=0.0 

QE-0  0 

DO  498  L-'l  ,  NSOIL 
498  CP  <  L ) =0  .  0 

497  APS-AE 

ITS=K 1+ 1 
XSS=0 . 0 
RETURN 
END 
C 

SUBROUTINE  Y ALIN ( UST , DFT , SNU , SUBW , S , G , VEL , HYR , CEE ) 

C 

C  THIS  SUBROUTINE  COMPUTES  THE  TRANSPORT  RATE  OF  NON  COHESIVE 
C  SEDIMENTS  USING  THE  BEDLOAD  FORMULA  DEVELOPED  BY  YALIN<1963> 

C 

CALL  SHIELD < UST , DFT , SNU, SUBW, TC) 
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Y=UST*UST/< <S-i . ) *G*DFT ) 

YC=TC/(SUBU*DFT) 

IF  <  Y-YC  >  4,4,5 

4  CEE=0 . 0 
GO  TO  6 

5  SEGMA=Y/YC-1 . 0 

A=2 . 45*SQRT(YC)/S**0  .  4 

X1  =  0 . 63S*DFT#UST*S*SEGMA 

X2= 1 . -ALOG ( i . +A*SEGMA ) / <  A*SEGMA ) 

CEE=X1*X2/(VEL*HYR> 

6  CONTINUE 
RETURN 
END 

C 

SUBROUTINE  YANG  <  DFT  ,  UST  ,  SNU ,  V ,  SL.P  ,  W ,  S ,  CEE  ) 

C 

C  THIS  ROUTINE  COMPUTES  THE  TRANSPORT  RATE  OF  NONCOHESIVE  SEDIMENTS 
C  USING  THE  TOTAL  LOAD  FORMULA  DEVELOPED  BY  YANG< 1973) 

D=DFT 

A=UST*D/SNU 

IF  <  A . GE .70.0)  GO  TO  7 

VCU=2 . 5/ <  ALOG 1 0  <  A ) -0 . 06>  +  0 .66 

GO  TO  8 

VCW-2 . 05 

ESP=V*SLP/W-VCW*SLP 
IF  ( ESP )  9,9,10 
CEE  =  0 . 0 
GO  TO  11 

0  FI =5 . 4 35- 0  286*AL0G 1 0 ( U*D/SNU >  -0  . 457*ALDGi 0 ( UST/W > 

F2=l .799-0 . 4 0 9* ALOG 1 0 ( W*D/SNU )  -  0  . 3 14* ALOG 1 0 (UST/U ) 
F3~ALOG10(ESP) 

E=Fi+F2*F3 
C=10 . 0**E 

CEE=C*S/(S-1 .  )*1  . 0E-6 
1  CONTINUE 
RETURN 
END 

SUBROUTINE  L AUR SN ( DFT , UST , S , V , GAMA , SUBW , SNU , DM , Y , G , W , CEE > 

THIS  SUBROUTINE  USES  THE  TOTAL-LOAD  FORMULA  DEVELOPED  BY 
LAURSEN( 19S8) 

D  =  DFT 
DP  =  DM 

CALL  SHIELD < UST , DFT , SNU, SUBW, TC> 

UU=ALOG10 (UST/W) 

IF(UU.GE . 0 . 40)  GO  TO  1 

FF=1 . /(12 . 32-10 . S*EXP ( 0 . 047*<UU+2 . > ) ) 

GO  TO  4 

IF < UU . GT .1.5)  GO  TO  2 
FF=2.04S*UU+0.942 
GO  TO  4 

IF  <  UU , GE .2.2)  GO  TO  3 

FF=3 . 38+SQRT ( 1 . 416- ( UU-2 . 51 >**2  .  ) 

GO  TO  4 

FF=0 . 26*UU+3 . 953 
FL-- 1  0 . 0**FF 

TO= ( (GAMA* V*V ) / ( G*58 . ) ) * ( DP/Y ) **0 . 3333333 
IF(TO  LT.TC)  GO  TO  5 

CEE  =  0  01*(  (TO/TO-l  .  )*FL*(D/Y)**1  .  1666666 
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GO  TO  6 
CEE=0  .  0 
CONTINUE 
RETURN 
END 

SUBROUTINE  SHIELD < UST , DFT , SNU , SUBU , TC > 

THIS  SUBROUTINE  COMPUTES  THE  CRITICAL  BED  SHEAR  STRESS  DERIVED 
FROM  SHIELDS'  FUNCTION 
REY=UST#DFT/SNU 
IF<REY . GT . 10.0)  GO  TO  1 
TC  =  0  . 08*SUPW*DFT/REY**D . 4 
GO  TO  3 

IF (REY . GT . 500 . 0  >  GO  TO  2 
TC  =  0 . 022*SUBW*DFT*REY**0  .  it 
GO  TO  3 

TC=0 . 06*Sl)BU*DFT 

CONTINUE 

RETURN 

END 

SUBROUTINE  SETVEL < D , T , S , W > 

THIS  ROUTINE  COMPUTES  THE  SETTLING  VELOCITY  OF  SEDIMENT  PARTICLES  OF 
ANY  DENSITY.  THE  ROUTINF.  ASSUMES  A  SHAPE  FACTOR  OF  0  7  AND  INTERPOLATES 
THE  VALUES  TABULATED  BY  THE  SUBCOMMITTEE  ON  SEDIMENTATION,  INTER¬ 
AGENCY  COMMITTEE  ON  WATER  RESOURCES,  1957. 

DIMENSION  A  ( 6 , 1 1 ) ,Z<2> 

DATA  A<l,l),A<l,2),A(l,3>,A(l,4>,A(l,S>,A<i,6>,A<i,7>,A<l,8>, 

$A<1,9> ,A(1,10) ,A<i,il)/ 

$0.04,0. 06,010, 0.20, 0.40, 0.80, 1.50,200, 3. 00, 700, 10. 00/ 

DATA  A<2,1>,A<2,2),A<2,3>,A<2,4>,A<2,5>,A<2,6>,A<2,7),A<2,8>, 

♦A (2, 9) ,A<2,10) , A<2 , 1 1 )  / 

$0 . 10,0 .24 ,0 . 60 ,1 .80 ,4 . 60 ,9 .50 , 16. 10 ,19. 90 ,25 .30 ,39 .50 ,44  00/ 

DATA  A(3,1),A<3,2),A(3,3),A<3,4>,A(3,5),A(3,6),A(3,7>,A(3,8), 

$A  <  3 , 9  > ,A(3, 10) ,A<3,11)/ 

$0 . 14,0 .32,0  76,2 .20 ,5 . 30 ,10 .50 , 16.90 ,20 .30 ,25 .60 ,39  50 ,44 .00/ 

DATA  A<4,1),A<4,2>,A(4,3),A<4,4),A<4,S>,A(4,6),A<4,7),A(4,8), 

$A  (4 ,9)  ,  A  <  4 , 1. 0  >  ,  A  ( 4 , 1 1 )  / 

$0 . 18,0 . 40 ,0 .92,2 .50 ,5 . 80 ,11 . 00 ,17 .SO ,20 .70 ,25.90 ,39  SO ,44 . 00/ 

DATA  A( 5 , 1 ) ,  A  ( 5 , 2 ) ,A(5,3) ,A(S,4) ,A(5,S> ,A<5,6) , A< 5 , 7 ) , A< 5 , 8 ) , 

$A<5 ,9) ,  A  ( 5 , 1 0  ) , A < 5 , 1 1 > / 

$0 .23,0 .49,1 . 10 ,2 .85,6.30 ,11 .60 ,17.90,21 . 10 ,26.20,39.50 ,44. 00/ 

DATA  A<6, 1 ) ,A<6,2) ,  A  <  6 , 3  ) ,  A  <  6 , 4  > ,A<6,5> , A < 6 , 6 ) ,A<6,7) ,A(6,8) , 

$A<6,9) ,A(6 , 10) ,A<6,11 >/ 

$0 .29,0 .57, 1 .26,3.20 ,6.70 ,12. 0  0,18  10,21 .50,26  SO ,39 .50 ,44 .00/ 
VSC(T)=1 . 41E-5  -3.48E-7«(T-10 . )  +5 . 0 OE-9* ( T-l 0 . ) *< T- 15 . ) 

$  +2.67E-11*<T-10 . >*<T  — 15  . >*(T-20 . >  -4  .  OOE-12*(T-10 . >* 

$  (T-15. )*<T-20. )*<T-25. > 

IF <D . GE . 0 .040)  GO  TO  2 
VISC=VSC(T) 

SS  =  0 .  0  009669*D*D/VISC 
GO  TO  18 

INTERPOLATION 
CONTINUE 
Q=T/10 . 
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KT=Q+i . 

PT=Q-KT+1 . 

DL=ALOG10(D) 

DO  10  J=i,10 

IF (D. LE . A( 1 , J> )  GO  TO  20 
10  CONTINUE 

20  J=J-i 

C=ALQGi 0 ( A  <  1 ,  J  )  ) 

E=ALOG10<A<l,J+i>> 

PD=><DL-C)/<E-C> 

DO  50  L=i,2 
I=L+KT 

50  Z<L )=< i . -PD)*ALOGl 0 <  A  <  I , J >  > +PD*ALOGi 0 ( A ( I , J  +  i) > 
R=<1  -PT)*Z<i)+PT*Z<2> 

SS—10  HR 

C  ADJUSTING  SETTLING  VELOCITY  FOR  SPECIFIC  GRAVITY 

18  W=SS*SQRT< (S-l ,0)/l .65)/30.5 

RETURN 
END 
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